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ABSTRACT 


Analytical  studies  combining  simple  one~dimensional  models  with  multidimensional 
computer  codes  have  been  performed  to  investigate  charging  characteristics  of  an 
emitting  probe  passing  through  a  partially  ionized  neutral  plasma.  The  effort  is  directed 
at  predicting  probe  potentials  as  a  function  of  emitted  electron  and  ion  currents  for  a 
variety  of  ambient  plasma  and  neutral  densities  and  temperatures  corresponding  to  those 
found  at  altitudes  between  150  and  400  km,  where  recent  sounding  rocket  experiments  have 
been  performed.  Calculations  have  been  made  for  the  early-time  transient  potentials  and 
probe  surface  currents  as  well  as  for  late-time  or  steady-state  potential-current  (l-V) 
probe  characteristics.  Efforts  have  been  directed  at  explaining  the  apparently 
anomalously  low  probe  potentials  and  local  maxima  in  probe  l-V  characteristics  which 
appear  to  occur  only  under  conditions  of  electron  emission. 

Various  possible  effects  previously  identified  were  investigated,  such  as  space-charge 
limiting  of  the  emitted  beam,  enhanced  local  plasma  produced  by  the  emitting  beams, 
rocket  probe  velocity,  the  earth's  magnetic  field,  and  sheath  dynamics  involving  both 
space-charge  limiting  and  orbit  limiting  of  the  return  currents  to  the  probe  which  act  to 
neutralize  the  beam  current.  Most  of  these  effects  were  investigated  in  the  early-time 
regime  (t  <  15  us)  with  the  two-dimensional  ABORC  computer  code,  and  virtually  all  of 
these  effects  except  the  enhanced  local  plasma  tend  to  increase  the  probe  potential  for  a 
given  beam  current.  None  appears  to  account  for  a  local  maxima  in  potential. 

During  the  course  of  this  effort,  ionization  of  the  neutral  gas  by  the  returning  probe 
current  (as  opposed  to  the  escaping  beam  current)  was  postulated  as  an  important 
mechanism  in  the  process.  Analytical  and  numerical  calculations  indicate  that  this 
mechanism  can  be  important  in  producing  anomalously  low  potentials  for  electron  current 
emission  but  not  for  ion  current,  apparently  in  agreement  with  rocket  observations.  It  also 
can  account  for  a  local  maximum  in  l-V  characteristics.  Calculations  with  the  ABORC 
code  in  the  early-time  regime  indicated  only  a  small  effect  due  to  ionization,  although 
potentials  had  reached  a  peak  and  were  beginning  to  decline.  Since  late-time  calculations 
could  not  be  performed  in  a  practical  manner  with  the  ABORC  code,  a  one-dimensional 
analytical  model  was  formulated  to  include  effects  of  both  space-charge  limiting  and 
ionization.  The  simplified  steady-state  model  showed  a  considerable  reduction  in 


predicted  potentials  for  electron-beam  emission  with  ionization,  demonstrating  that 
ionization  of  the  background  gas  is  a  major  factor  in  determining  the  probe  potential  in 
regimes  where  a  significant  gas  density  exists.  For  conditions  that  were  investigated  in 
some  detail  (Te  *  1000° K,  Np|asma  =  103/cm2,  Nneutra|  =  5  x  lO^/cm3),  the  probe 
potential  initially  increased  with  increasing  emission  current.  As  the  rate  of  ionization 
became  significant  with  increasing  emission  current,  the  probe  potential  reached  a 
maximum  and  then  decreased  at  higher  emission  currents.  For  high  neutral  densities 
(~10  /cm  )  and  relatively  high  emission  currents  ( — 0 .1  A),  the  probe  potential  appears  to 
become  pinned  to  the  assumed  ionization  threshold  (~40  to  50  eV).  It  should  also  be  noted 
that  the  potential  curve  with  respect  to  spatial  position  is  very  flat  near  the  probe,  so  that 
the  potential  difference  between  the  probe  and  a  point  located  a  few  meters  or  less  from 
the  probe  would  be  much  less  than  the  probe  potential  with  respect  to  infinity. 


1.  INTRODUCTION 


1.1  BACKGROUND 

When  an  emitting  probe  is  inserted  into  a  partially  ionized  neutral  plasma,  it  will 
attain  a  potential  relative  to  the  ambient  plasma  that  is  a  function  of  several  parameters, 
among  which  are  the  size  and  shape  of  the  probe,  probe  velocity,  plasma  temperature  and 
density,  neutral  gas  density  and  species,  magnetic  field,  beam  current  and  voltage,  and 
beam  polarity-  The  relative  importance  of  these  parameters,  the  functional  dependence  of 
the  probe  current-voltage  (l-V)  characteristics  on  these  parameters,  and  the  way  in  which 
the  beam  couples  to  the  ionosphere  are  not  yet  well  understood. 

To  investigate  these  effects,  a  number  of  experiments  involving  sounding  rockets  in 
the  upper  ionosphere  have  been  performed  in  which  the  rocket  potentials  have  been 
measured  for  various  ion  and  electron  beam  conditions.  The  literature  describing  these 
experiments  is  included  in  a  companion  literature  summary  report  (Kef-  1).  Kesults  of 
these  experiments  indicate  relatively  high  probe  potentials  (~85U  V)  for  relatively  low  ion 
beam  currents  ( —  1 0  ^  A),  and  relatively  low  probe  potentials  (~100  V)  for  relatively  high 
electron  currents  (~0.1  A).  In  addition,  there  is  some  evidence  that  a  local  maximum  in 
probe  potential  as  a  function  of  emitted  electron  current  has  been  observed  (Kef.  2). 

Numerous  analytical  and  computer  calculations  have  been  made  in  an  attempt  to 
explain  the  observed  variation  in  probe  potentials  with  altitude,  emission  current  and 
species,  and  rocket  orientation  with  respect  to  the  geomagnetic  field.  Much  of  this 
literature  is  also  described  in  Keference  1. 

It  appears  that  the  probe  potential  for  ion  emission  can  be  explained  on  the  basis  of 
Langmuir  probe  theory,  including  space-charge  and/or  orbit  limiting,  the  theories  of  which 
are  described  in  numerous  references  (see  Kefs.  3-5).  These  theories  have  to  be  somewhat 
modified  to  account  for  the  geomagnetic  field  (Kefs.  6-7). 

The  simple  Langmuir  probe  theory  does  not  account  for  the  low  observed  potentials 
associated  with  electron  emission,  however.  Most  of  the  postulated  effects  which  have 
not  been  treated  in  detail,  such  as  the  geomagnetic  field,  rocket  velocity,  and 
mult idimensionalitv,  would  appear  to  result  in  higher  probe  potentials  rather  than  in  lower 
potentials.  O'Neil  et  al.  (Ref.  6),  using  Langmuir  probe  theory  with  an  enhanced  steady- 
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state  electron  concent  rat  ion  near  the  probe  produced  by  ionization  of  the  ambient 
background  by  the  emitting  beam,  was  able  to  arrive  at  relatively  low  potentials  for 
electron  beam  emission. 

This  explanation  is  not  entirely  satisfying  for  several  reasons,  particularly  for  high- 
energy  emitted  beams,  First,  any  ionization  produced  by  the  beam  would  be  highly 
localized  around  the  beam,  and  the  one-dimensional  Langmuir  theory  used  in  Reference  8 
assumes  a  uniform  enhanced  distribution  around  the  entire  probe.  Second,  and  just  as 
important,  the  return  current  to  the  probe  equals  the  emitted  current  in  steady  state  and 
is  of  lower  energy,  often  a  few  hundred  elect ronvolts.  The  ionization  rate  produced  by 
these  low-energy  returning  electrons  would  be  substantially  higher  than  that  occurring 
from  the  high-energy  emitted  electrons.  Thus,  ionization  by  the  return  current  would 
dominate  the  ionization  produced  by  the  emitted  current. 

It  is  expected  that  the  probe  response  would  be  considerably  different,  depending  on 
which  mechanism  dominates.  If  ionization  by  the  beam  dominates,  the  return  current  to 
the  probe  would  be  largely  confined  to  the  region  about  the  beam  and  the  net  beam 
current  leaving  the  probe  would  be  significantly  smaller  than  the  emitted  beam  current. 
Currents  flowing  on  the  rocket  probe  will  be  correspondingly  small.  If  the  return  current 
is  relatively  uniform  through  the  plasma,  the  net  beam  current  will  be  approximately  equal 
to  the  emitted  beam  current,  and  substantially  higher  replacement  currents  will  flow  on 
the  rocket  probe. 

1.2  APPROACH 

The  main  thrust  of  the  present  effort  is  to  explore  the  mechanisms  involved  in  beam 
coupling  to  the  ionosphere,  to  quantify  the  probe  l-V  characteristics  as  a  function  of 
relevant  physical  parameters,  and  to  identify  a  plausible  mechanism  which  might  explain  a 
local  maximum  in  the  l-V  characteristics  for  an  electron-emitting  probe. 

During  the  course  of  these  efforts,  it  became  apparent  that  the  ionization  of  the 
neutral  background  gas  by  electrons  returning  to  the  vehicle  to  balance  the  emitted 
current  would  play  a  dominant  role  in  determining  the  probe  potential.  The  main  question 
is  whether  this  process  can  produce  the  required  magnitude  of  reduction  in  probe 
potentials  at  the  densities  of  neutral  plasmas  corresponding  to  existing  data,  and  whether 
the  functional  dependence  of  the  potentials  on  the  various  physical  parameters  can  be 
explained  by  this  mechanism. 
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A  two-fold  approach  was  adopted  with  the  aim  of  providing  some  quick  order-of- 
magnitude  estimates  of  the  effects,  along  with  a  longer-term  effort  that  would  quantify 
the  effects  in  more  detail,  perhaps  in  multidimensional  geometries  if  necessary.  These 
approaches  included  simple  one-dimensional  analytical  steady-state  models  using  Poisson 
equations  with  ionization  of  the  neutral  background  by  returning  electrons,  and  a  two- 
dimensional  time-dependent  computer  code  with  space-charge  limiting,  ionization,  and 
magnetic  field  effects.  Near  the  end  of  this  effort,  a  quasi-static  two-fluid  model  of  the 
electron  and  ion  gas  was  employed  in  an  attempt  to  overcome  computational  difficulties 
of  the  particle-pushing  computer  code  and  obtain  two-dimensional  time-dependent 
information. 

The  principal  analytic  treatment  consisted  of  a  steady-state  solution  of  the  Langmuir 
space-charge  problem  with  ionization  in  spherical  symmetry.  Because  of  the  spherical 
symmetry,  these  solutions  cannot  rigorously  consider  the  effects  of  the  geomagnetic  field, 
the  velocity  of  the  probe,  or  the  dynamics  of  the  emission  beam  itself.  In  Reference  3, 
Lam  presented  a  solution  to  the  problem  of  an  emitting  spherically  symmetric  probe 
without  ionization.  This  is  one  of  the  simplified  models  that  predicted  too  large  potentials 
for  electron  emission.  To  gain  greater  insight  into  the  factors  that  influence  the 
magnitude  of  the  probe  potential,  this  spherically  symmetric  model  was  reexamined  under 
this  program,  with  a  slightly  different  outer  boundary  condition  than  the  one  used  by 
Lam.  A  series  solution  was  found  for  the  potential  in  the  sheath  depletion  region  which 
agreed  very  well  with  Lam's  numerical  integration  (Section  3).  Unfortunately,  it  does  not 
appear  practical  to  include  ionization  effects  and  produce  purely  analytical  solutions. 
Consequently,  the  equations  were  solved  numerically  for  the  sheath-region  potential  in  a 
form  that  could  later  be  extended  to  include  ionization.  As  would  be  expected,  the 

numerical  results  from  the  computer  code  without  ionization  are  in  reasonable  agreement 

with  the  analytic  results  of  Lam  and,  therefore,  with  the  series  solution.  Consequently, 
these  calculations  also  predict  too  large  potentials  for  electron  emission  without 
ionization. 

When  ionization  is  added  to  the  problem,  the  situation  becomes  much  more  com¬ 
plicated  because  the  charge  density  at  a  particular  radius  from  the  probe  is  not  just  a 
function  of  the  potential  at  that  position,  as  it  is  in  the  simple  models  without  ionization, 
but  is  a  function  of  the  ionization  rate  at  all  other  points  in  space  and  of  the  potential 
difference  from  the  point  of  interest  to  those  positions.  In  other  words,  the  charge 

carriers  at  a  given  location  retain  some  memory  of  where  and  at  what  potential  they  were 

created. 
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To  solve  this  problem,  the  previously  developed  numerical  technique  was  extended  to 
include  ioni/ation,  as  described  in  Section  3.4.2.  The  solution  to  the  problem  with 
ionization  is  computationally  quite  sensitive,  especially  for  mean  free  ionization  distances 
less  than  10  to  100  m  and  large  emission  currents.  The  result  is  that  the  relatively  simple 
iterative  procedure  used  in  the  code  would  not  always  reach  a  fully  converged  solution. 
The  difficulty  undoubtedly  could  be  improved  with  more  refined  iteration  algorithms. 
However,  at  worst,  the  calculated  results  tend  to  oscillate  about  the  apparently  correct 
solution.  The  results  clearly  indicate  that  even  a  relatively  low  background  gas  density 
(relatively  large  ionization  distance)  can  greatly  reduce  the  probe  potentials  from  the 
predicted  values  without  ionization. 

An  alternate  approach  attempted  late  in  this  effort  was  a  quasi-static  time-dependent 
two-fluid  computer  model  of  the  electron  and  ion  gas,  programmed  in  cylindrical ly 
symmetric  coordinates  using  finite-difference  techniques;  this  code  was  a  relatively  small 
modification  to  a  previously  existing  code.  It  was  to  be  utilized  in  an  attempt  to 
overcome  the  computational  difficulties  which  arise  in  all  particle-pushing  codes,  such  as 
ABORC,  when  the  particle  densities  become  too  large.  One  of  its  advantages  over  the 
steady-state  models  is  that  it  could  give  the  time  evolution  of  the  probe  potential.  In 
addition,  it  can  treat  more  realistic  cylindrical  probes.  Some  early  results  gave  promise 
that  this  code  could  be  useful,  especially  in  the  high-density  regimes.  Unfortunately,  no 
specific  results  were  obtained  at  the  time  this  report  was  prepared.  Such  an  approach 
appears  promising  and  will  be  pursued  further.  The  physics  and  status  of  this  technique 
are  documented  in  Appendix  A. 

The  computer  code  effort  was  initiated  by  a  review  of  the  known  computer  codes 
thought  to  be  potentially  useful  for  space  plasma  problems  in  general  and  the  present 
problem  in  particular  (Ref.  1).  The  "two-and-a-half "-dimens ion,  cylindrically  symmetric 
computer  code  ABORC  (Ref.  9)  contains  many  of  the  desirable  features  for  the  present 
problem  and  is  reasonably  economical  to  run.  This  code  was  modified  to  include  the 
earth's  magnetic  field,  an  initial  ionized  background  plasma,  and  ionization  of  the  neutral 
gas  by  the  emitted  electron  beam  and  the  replacement  plasma  electrons.  The  advantage 
of  this  code  is  that  its  physics  are  quite  rigorous,  with  a  minimum  of  simplifying 
assumptions  and  approximations.  Its  main  disadvantage  is  that,  for  computational  reasons, 
it  is  practical  to  use  this  code  only  for  early  times  (t  <  50  jjs)  after  the  emission  beam  is 
turned  on.  In  this  time  regime,  it  is  adequate  to  simulate  the  positive  ions  as  an  immobile 
uniform  background  charge  density,  although  that  restriction  is  not  inherently  necessary  in 
the  ABORC  code.  1  he  modifications  that  were  made  to  the  ABORC  code  specifically  for 
this  program  are  described  in  Section  4,  along  with  the  computations  and  results. 
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1*3  REPORT  ORGANIZATION 

The  initial  literature  survey  of  analytic  techniques,  numerical  computer  codes,  and 
rocket  orobe  experiments  is  contained  in  a  separate  companion  volume  (Kef.  1).  The 
ionospheric  environments,  including  neutral  and  ion  densities  and  temperatures  as  a 
function  of  altitude  and  time  of  day,  are  summarized  in  Section  2,  along  with  a  sample  of 
the  pertinent  rocket  probe  data  obtained  for  both  ion  and  electron  emission. 

The  steady-state  analytic  approach  using  Poisson's  equation  is  described  in  Section  3. 
Comparisons  to  early  work  by  Lam  (Ref.  3)  in  the  absence  of  ionization  are  also  included. 
The  principal  results  of  this  report,  describing  the  effect  of  ionization  of  the  neutral  gas 
by  returning  electrons  and  the  local  maxima  in  the  l-V  characteristics,  are  also  described. 

The  numerical  results  using  the  two-dimensional  time-dependent  ABORC  code  are 
described  in  Section  4.  Comparisons  are  made  with  the  one-dimensional  numerical  work  of 
Rothwell  (Ref*  10).  The  effects  of  ionization  on  the  early-time  response  are  reported,  as 
well  as  the  effects  of  the  magnetic  fields.  One-dimensional  and  two-dimensional  results 
are  also  compared. 

The  main  results  of  the  effort  are  summarized  in  Section  5.  A  short  description  of  an 
alternate  numerical  approach  thought  to  be  capable  of  describing  the  time-dependent 
behavior  of  the  rocket  probe  in  two  dimensions  is  summarized  in  Appendix  A. 
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2,  ENVIRONMENTS  AND  DATA  FOR  IONOSPHERIC  PROBE  EXPERIMENTS 


This  section  briefly  describes  the  ionospheric  parameter  ranges  of  interest  to  the 
present  work,  and  surveys  experimental  data  obtained  from  charged  particle  beam 
excitations  by  probes  at  various  altitudes.  Only  the  parameters  relevant  to  the 
calculat iona I  modeling  for  probe  potential  are  considered.  Additional  details  of  the 

environment  are  discussed  in  Reference  1. 

At  high  altitudes,  the  sun's  radiation  causes  appreciable  photoionization  of  the 

atmosphere.  At  low  gas  pressures,  recombination  of  the  electrons  and  ions  is  slow  enough 

that  high  electron  concentrations  can  exist  even  through  the  night.  Figure  1  shows  typical 

day  and  night  electron  densities  as  a  function  of  altitude  for  the  extremes  of  the  sunspot 

cycle.  Above  the  F2  region,  the  electron  density  monotonically  decreases  out  to  several 

earth  radii.  Beyond  several  earth  radii,  at  the  outer  edge  of  the  protonosphere,  the 

electron  densities  are  determined  by  solar  wind  or  interplanetary  plasma.  For  100-  to  500- 

3  6  T 

km  ranges,  the  electron  density  varies  from  10  to  10  e/cm  .  The  electron  temperature 
can  be  as  large  as  1  eV,  as  discussed  in  Reference  1.  Other  less  relevant  plasma 

characteristics  for  the  present  purposes  are  also  described  there. 

The  density  of  neutral  species  is  a  strong  function  of  height  and  exospheric 

temperature.  Figure  2  shows  different  species  concentrations  as  a  function  of  altitude  for 

three  different  temperatures.  The  "geometric  height"  corresponds  to  the  local  altitude  of 

the  distributions  shown  for  a  given  latitude,  and  the  "geopotential  height"  corresponds  to 

the  altitude  at  a  given  latitude  of  the  same  isopotential  surface  (Ref.  11).  A  reasonable 

13  3 

maximum  value  is  approximately  10  neutrals/cm  when  all  the  species  are  added 
together  at  the  100-km  height.  It  is  reasonable  to  lump  the  species  together  in  the  present 
computer  model  because  the  cross  sections  for  ionization  by  electrons  are  similar  in 
magnitude  and  their  dependence  on  electron  energy.  Notice  also  that  a  mean  molecular 
mass  of  30  for  the  neutrals  gives  the  ions  a  mass  of  approximately  55,000  times  the 
electron  mass  (corresponding  to  a  velocity  ~200  times  slower  for  a  given  energy).  Also, 
variations  as  a  function  of  latitude  or  isopotential  surfaces  are  not  of  prime  concern  here, 
but  rather  scoping  the  parameter  range  is  important. 
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ELECTRON  CONCENTRATION  (e/cm3) 


Figure  1.  Example  of  normal  electron  distributions  at  the  extremes 
of  the  sunspot  cycle  at  geomagnetic  latitudes 
of  30  to  40*.  Curves  are  from  Reference  12. 


Ionospheric  probe  data  from  several  experimenters  is  summarized  in  Table  1  and 
Figure  3,  The  table  lists  the  experimental  and  environmental  parameters.  Both  electron 
and  positive  ion  beams  of  various  currents  (I)  and  energies  (<?)  were  employed.  A 
considerable  altitude  variation  (h)  resulted  in  a  wide  range  of  ambient  electron  densities 
(n).  The  temperatures  and  densities  of  the  background  environment  are  not  provided  in  the 
references  in  several  cases.  Probable  values  from  environmental  graphs  presented 
previously  are  indicated  by  a  question  mark  in  the  table  for  those  experiments.  Perhaps 
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NUMBER  DENSITY  (m*3)  NUMBER  DENSITY  (m‘ 


GEOMETRIC  HEIGHT  (km)  GEOMETRIC  HEIGHT  (km) 


0  500  1000  1500  2000  O  500  1000  1500  2000 

GEOPOTENTIAL  HEIGHT  (km)  GEOPOTENTIAL  HEIGHT  (km) 


GEOMETRIC  HEIGHT  (km) 


Figure  2.  Number  densities  of  individual  atmospheric  constituents  as  a 

function  of  height  for  three  representative  exospheric  temperatures. 
The  mean  molecular  mass  as  a  function  of  height  is  shown  for 
various  exospheric  temperatures  in  the  lower  right  diagrams. 

Curves  are  from  Reference  11. 
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the  most  striking  result  of  the  data  ts  the  tendency  of  the  probe  to  linger  at  100  V  or  less 
(as  measured  to  ~3  m  out  from  the  probe  in  some  cases)  for  electron  beam  emission.  Ion 
emission  at  a  much  lower  beam  current  than  any  of  the  electron-beam  cases  resulted  in  a 
potential  at  least  10  times  higher  for  similar  environmental  parameters.  Possible 
explanations  for  this  behavior  have  included  the  background  particle  mass  differences,  the 
geomagnetic  field,  and  the  electron-neutral  ionization  cross  section  (which  is  maximum 
near  100  eV).  The  last  item  appears  to  be  the  major  effect  limiting  the  probe  potential,  as 
reported  in  Section  3.  The  ions  have  a  much  lower  ionization  cross  section  than  the 
electrons,  and  allow  the  potential  to  increase  further  than  an  electron-beam  case.  The 
analytic  theory  of  Lam  (Ref.  3),  which  does  not  include  ionization  effects,  was  evaluated 
and  found  to  be  in  reasonable  agreement  with  the  data  for  the  ion  case  (see  the  data  points 
in  Figure  3).  The  analytical  results  are  discussed  further  in  Section  3. 


Table  1.  Ionospheric  Probe  Experimental  Data  Summary 


Experimental  Parameters 

Environment  Spec. 

E  xperimenter 

p.W., 

(A) 

(eV) 

fy* 

h 

(km) 

n 

(#/cm3) 

T 

(eV) 

- 

Rocket/Size 

Date 

Reference 

Cohen 

e” 

several 

3x103 

100 

120 

103? 

0.1? 

? 

1979? 

Private  comm, 
to  H»  Linnerud 

Cohen 

«’ 

0.1 

45x103 

100 

350 

2*105? 

0.1? 

? 

1979? 

3/80 

R,  O'Neil 

e 

0.6 

2.5  x103 

3-30 

80-120 

103? 

0.14 

Precede 

6  x  0.23  m 

10/74 

ICR  7/1/78 
(Re?.  8) 

Arnoldy/Winckler 

e” 

0.07 

4x10* 

3-5 

350? 

2x105 

0.1? 

Echo  ill 

>1978 

Arnotdy  8r  Winck- 
ler,  UNH  rpt  78 
(Ref.  13) 

Echo  IV,  Ref.  14 

Cohen 

e~ 

0.01 

90 

90 

150? 

0.1? 

White  Sands 
3  x  0.38  m 

1/21/78 

GRL  6/79,  Cohen 
(Ref.  15) 

Cohen 

Xj(Z«54) 

1.2x10  ^ 

2x103 

850 

100-400 

103 

0.1? 

White  Sands 
3  x  0.38  m 

1/21/78 

GRL  6/79,  Cohen 
(Ref.  15) 

f  rench- Russian 

e"? 

0.5 

1.5x10* 

100 

? 

? 

? 

ARAKS 

1974 

Gendrin 
(Ref.  20) 

RE-03225 


ELECTRON  DENSITY  (e/cm3) 


Figure  3.  Experimental  data  of  Cohen  (Ref*  15)  for  ion  beam  emission 
from  a  rocket  probe 


3.  SPHERICALLY  SY  M  METRIC  STEADY-STATE  MODEL 


3.1  LAWS  MODEL 

In  Reference  3,  Lam  presented  a  steady-state  analysis  for  a  spherically  symmetric 
emitting  probe  without  ionization*  His  model  is  an  extension  of  the  l  angmuir-Chiids 
solution,  using  a  different  outer  boundary  condition.  He  divided  the  radial  distance  from 
the  probe  into  three  regions:  (1)  a  quasi-charge-neutral  outer  region  where  the  electron 
and  ion  densities  are  approximately  equal,  (2)  an  inner  sheath  region  where  the  charge 
species  that  is  repelled  from  the  probe  has  essentially  zero  density,  and  (3)  an 
intermediate  transition  region. 

It  is  useful  to  repeat  his  mathematics  for  the  sheath  region  since  it  is  identical  to  the 
L angmuir-Chiids  model  used  for  the  series  solution  (Section  3.2)  and  nearly  identical  to 
what  is  used  in  the  computer  model  without  ionization  in  Section  3.4.1. 

The  potential  <\>  at  any  point  is  governed  by  Poisson’s  equation: 


v2 


♦ 


-P/En 


(3-1) 


where  cQ is  the  permittivity  of  free  space  and  pis  the  charge  density  of  the  attacted 
species,  since  the  repelled  species  is  assumed  to  be  completely  depleted  from  the  sheath 
region.  The  density  of  the  emitted  particles  is  ignored  here  since  their  velocities  are 
usually  much  greater  than  the  returning  particles  and,  hence,  their  densities  are  much  less. 

In  steady  state,  the  returning  current  must  be  continuous  and  equal  to  the  emitted 
current: 


I  =  - pvA  =  -4 Trpvr2  ,  (3-2) 

where  positive  I  is  an  outward  going  current,  v  is  the  absolute  velocity  of  the  returning 
particles,  and  A  is  the  area  of  an  imaginary  shell  located  at  r. 

Combining  Eqs.  3-2  and  3-1, 


In  the  sheath  region,  the  potential  will  often  be  considerably  greater  than  the  initial 
thermal  energy  of  the  charged  particles.  In  that  case,  the  particle  velocity  v  can  be 
obtained  by  equating  the  kinetic  energy  of  the  particle  to  the  electrostatic  potential 


energy: 


j  m  v2  =  e|  4-1  ,  (3-4) 

where  m  is  the  mass  of  the  particle  and  e  is  the  magnitude  of  the  electronic  charge.  This 
assumption  is  sometimes  called  the  'free-fall*  approximation  for  the  attracted  species. 
Combining  Eqs.  3-4  and  3-3  gives 


4i,eo  f2  4*1 

In  the  sheath  region,  Lam  defines  a  dimensionless  potential, 


and  a  dimensionless  inverse  distance. 


(3-5) 


(3-6) 


T 


ro/r  ' 


(3-7) 


where  rQ  is  the  outer  radius  of  the  sheath  region.  We  have  added  the  absolute  value 
symbols  in  Eq.  3-6  to  make  the  results  apply  to  either  electron  or  ion  emission.  In  terms 
of  these  variables,  Eq.  3-5  becomes 


2  92Jl  _  1_ 

8t2  /? 


(3-8) 


Lam  integrates  Eq.  3  -tt  numerically  across  the  depletion  region  using  a  boundary 
condition  at  r  =  r^  which  he  obtained  by  matching  numerical  solutions  of  his  differential 
equations  in  the  other  two  regions.  His  solution  is  reproduced  here  in  Table  2. 
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Table  2.  Tabulation  of  Function  F(t) 


T 

F(t)  (from  Lem) 

F(T)  (from  series  solution) 

1.0 

0.0000 

0.0000 

1.01 

0.00370 

0.00369 

1.02 

0.00929 

0.00927 

1.03 

0.01590 

0.01591 

1.04 

0.02327 

0.02322 

1.05 

0.03125 

0.03123 

1.10 

0.07775 

0.07767 

1.15 

0.13189 

0.1319 

1.20 

0.19131 

0.1912 

1.25 

0.25473 

0.2545 

1.295 

0.31456 

0.3144 

1.30 

0.32135 

0.3212 

1.S0 

0.610 

0.6102 

2.00 

1.41 

1.417 

3.00 

3.19 

3.173 

10.00 

16.4 

16.28 

100.0 

188.0 

188.17 

3.2  SERIES  SOLUTION  IN  SHEATH  REGION 

If  the  outer  two  regions  in  Lam's  analysis  are  ignored  and  <J>  is  taken  to  be  zero  at  r  = 
Tq  (as  is  done  in  the  computer  model  in  Section  3.4),  a  solution  to  tq.  3-8  can  be  obtained 
exactly  in  terms  of  a  series  solution: 


F  =  A.  Jin 


4/3 


+  A2 


in 


7/3 


t  ♦  in 


10/3 


T  ♦ 


(3-9) 


This  series  automatically  satisfies  the  two  boundary  conditions  that  the  potential  (F)  and 
the  electric  field  (3F/3r)  are  zero  at  the  outer  boundary  of  the  sheath  (t  =  1.0).  The 
constants  Aj  can  be  evaluated  by  substituting  tq.  3-9  into  Eq.  3-8  and  equating  equal 
powers  of  Jin  x  on  both  sides  of  the  equation.  The  result  is 


F=  1.7171  Jtn4/3  t  +  0.68684  Sjn/l  x  ♦  0.20605  £n10/3  i  *  0.04766  £n13/3  i 


(3-10) 


+  0.00894  fcn16/3  x  ♦  0.00141  £n19/3  x  ♦ 


The  series  in  Eq.  3-10  may  be  convergent  for  all  values  of  in  x  since  the  coefficients  of 
F/1.7171  in*^  x  shown  in  Eq.  3-10  are  smaller  than  those  in  the  expansion  of  (ev-1)/v,  v  - 


However,  we  have  not  rigorously  proven  that  the  higher  order  terms  are  convergent. 
C onsequent ly ,  this  series  will  only  be  used  for  0  <  in  t<  1,  where  it  is  clearly  convergent. 

If  the  ratio  of  the  sheath  radius  to  the  probe  radius  is  not  too  large  (that  is, 
r^K  <2. 718),  tq.  3-10  can  be  used  for  the  potential  profile  throughout  the  sheath  region. 
For  values  of  r^/K  >  2.718,  another  series  solution,  which  also  satisfies  Eq.  3-8,  can  be 
used: 


c  i,  ^  -1/2  t>  -3/2  „  -2  u  -5/2 

F  =1^  i  t  +  B4  t  +  B  t  +  T 


0  -3  M  -7/2  u  -4 

+  B7  T  +  Bft  T  ♦  t  +  .  .  . 


(3-11) 


Note  that  the  terms  and  t  ^  do  not  appear  in  this  series.  If  this  series  and  the  first 

derivative  of  F  with  respect  to  r  are  convergent  for  T  =  1.0,  Eq.  3-11  can  be  used 
throughout  the  sheath  region  from  r  =  r^  to  r  =  K.  The  coefficients  B.  (i  >3)  can  be 
obtained  in  terms  of  B^  and  l*2  by  inserting  Eq.  3-11  into  Eq.  3-8  and  equating  coefficients 
of  equal  powers  of  t  on  the  two  sides  of  the  equation.  and  B^  can  then  be  obtained 
using  the  boundary  conditions  for  the  potential  and  the  electric  field  at  the  outer  edge  of 
the  sheath  region  (both  =  0).  However,  if  these  series  are  not  sufficiently  convergent  at  T 
=  1.0,  Eq.  3-10  can  be  used  for  the  potential  profile  for  1  <  x<  2.718,  and  Eq.  3-11  can  be 
used  for  t  >  2.718.  I  he  coefficients  B^  and  B2  can  then  be  obtained  by  matching  the 
magnitude  and  first  derivative  of  Eq.  3-11  to  the  magnitude  and  first  derivative  of  Eq.  3- 
10  at  some  arbitrary  value  of  t  such  as  T  =  2.718,  where  Eq.  3-10  is  still  valid  and  Eq.  3-11 
is  more  convergent  than  at  x  =  1.0.  This  is  the  procedure  that  has  been  used  here.  The 
resulting  values  for  the  first  few  B  coefficients  and  the  analytic  relations  of  B-  (i  >3)  to  B^ 
and  B2  are  given  below. 

B1  =  1.9123 

B2  =  -3.1532 
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=  -U.0424 


B 

B 


7-'  h  B4 
#  =  ‘  63  B5 


«.  I.  H  u  b-5/2 
16  H2  B3  B1 


*  h  “3  's5/2 


ii  A 


0.0727 


The  value  of  F  vs-  T  from  these  two  series  solutions  (Eq.  3-10  for  1  <  2.718  and  tq.  3-11  for 
x>  2.718)  are  compared  to  Lam's  results  in  Table  2.  The  agreement  is  very  good  over  the 
whole  range  of  t  in  the  table. 


33  ANALYTICAL  RESULTS 


The  above  results  have  been  used  to  generate  some  useful  parametric  curves.  From 
Eq.  3-6,  the  probe  potential  (at  r  =  R),  relative  to  infinity,  is  given  by 


=  U<R> 


|l|2/3  F(r0/K)  . 


(3-12) 


For  relatively  large  emission  currents  I  (and  therefore  large  probe  potentials  <J>),  the  outer 
boundary  of  the  sheath  region  Tq  can  be  assumed,  with  acceptable  accuracy,  to  be  the 
radial  position  where  the  inward  thermal  current  of  the  attracted  carriers  just  equals  the 
emission  current: 

l  =  )A=-n0q  vth4n  ,  (3-13) 

where  vt^  =  (kT/2ir  m)1^2,  nQ  is  the  density  of  the  attracted  species  in  the  bulk,  q  is  its 
charge  (+  for  ions,  -  for  electrons),  kT  is  its  thermal  energy,  and  v  ^  is  the  average 
velocity  of  the  total  density  nQ  crossing  a  plane  in  one  direction.  If  one  considers  only  the 
half  of  n^  whose  velocities  are  directed  toward  a  given  plane,  their  average  velocity 
toward  the  plane  is 

/M  . 

tt  m 

Solving  Eq.  3-13  for  r^,  the  value  of  i  corresponding  to  the  probe  radius  R  is 


T 

P 


r„/R  = 


y  4  n  k' 


-i 


q  v 


th  0 


(3-14) 


Note  that  -l/q  is  always  a  [>ositive  quantity. 


I  or  purposes  of  numerical  evaluation,  we  take  the  typical  ion  mass  here  to  he  46,000 
tunes  the  electron  mass,  which  corresponds  to  a  mean  molecular  mass  of  about  25,  based 
on  figure  2  for  an  altitude  of  about  150  km.  With  this  mass,  the  potentials  become 


<)■  =  b.1  x  102  \L/S  \  ( r()/ 1< )  (3-15) 

and 

*=  -2.2  x  104  |l|2/3  F(r0/K )  (3-16) 


for  electron  anti  ion  emission,  respectively. 

It  is  convenient  to  specify  explicitly  the  function  F(tq/R)  in  terms  of  the  current  I. 
lo  do  this,  we  introduce  a  current  l(),  defined  as  the  current  due  to  the  thermal  electrons 
which  would  strike  the  probe  if  the  probe  potential  were  zero,  and  thus  the  bulk  plasma 
densities  and  velocities  extended  all  the  way  into  the  probe  and  radius  R: 


I 

o 


-4  ii  R  q  n 


0 


(3-17) 


I  aking  the  square  root  of  the  ratio  of  tqs.  3-13  and  3-17, 


o 


The  quantity  I  /K  is  shown  in  Figures  4  and  5  for  electrons  and  for  ions  with  a  mass 
m.  =  46,000  times  the  electron  mass  m  respectively,  as  a  function  of  n  for  various 
temperatures.  The  l-V  characteristics  for  the  probe  are  shown  in  Figures  6  and  7  for 
various  values  of  I  which  in  turn  depends  on  R,  T,  and  n. 

It  is  interesting  to  use  these  curves  to  predict  the  potentials  for  typical  experimental 
conditions  listed  in  fable  1  and  Figure  3.  For  electron  emission,  use  an  emission  current  I 
=  0.1  A  with  a  plasma  density  n  =  2  x  10^/cm^  and  a  temperature  1  =  0.1  eV,  and  assume 
an  effective  probe  radius  of  0.5  m.  From  Figure  4,  Iq/K  =  0.022  A/m  ,  so  I  =  5.5 
xIO  *  A.  from  Figure  6,  ^  -  900  V,  which  is  considerably  larger  than  the  potentials 

measured  by  Cohen  and  by  Arnoldy/  Winckler  for  comparable  conditions  (Table  1).  For  ion 

-  5  3  3 

emission,  use  I  =  1.2  x  10  A,  n  =  10  /cm  ,  l  =  0.1  eV,  and  R  =  0.5  m.  From  Figure  5, 
l()/R^  =  5  x  10  ^  A/m^,  so  l(f  =  1.25  x  10  ^  A.  From  Figure  7,  <p ^  -  300  V,  which  is  in  fair 
agreement  with  Cohen's  data  in  Table  1  for  ion  emission  (850  V),  considering  the 
uncertaint  ics  in  the  environment  conditions  and  the  average  value  for  R. 
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Figure  4.  Thermal  electron  current  vs*  electron  density  and  temperature 

For  the  theoretical  data  points  in  Figure  3,  an  effective  spherical  probe  radius  of 
0.56  m  was  used  to  give  the  same  total  probe  area  as  on  the  cylindrical  probe  used  in  the 
Cohen  experiments  (Ref.  15).  As  noted  later,  this  assumption  could  help  to  explain  why 
the  theoretical  values  are  somewhat  less  than  the  measured  values.  The  mass  of  the 
returning  ions  was  assumed  to  be  46,000  times  the  electron  mass,  consistent  with  Eq.  3-16. 
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Figure  5.  Thermal  ion  current  vs.  ion  density  and  temperature 


For  the  emission  current  of  8  p  A  in  Figure  3,  the  outer  radius  of  the  sheath  region  (Tq) 

is  4  m,  1.26  m,  and  0.4  m  for  bulk  ion  densities  of  10**,  10*,  and  loVcm3,  respectively. 

Equation  3-16  and  Table  2  (or  Figures  5  and  7)  then  yield  the  theoretical  potentials  from 

the  probe  to  infinity  shown  in  Figure  3  for  nQ  =  103  and  104/cm3.  It  should  be  noted  that 

3  3 

the  4  m  sheath  radius  for  10  /cm  is  outside  the  range  of  the  potential-measuring  device 
used  in  Reference  15.  Flence,  for  a  direct  comparison  to  the  experimental  data  in  Figure  3 
at  103/cm3,  the  calculated  potential  should  be  reduced  somewhat.  Since  the  value  of  r^  = 
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0.4  m  for  10  /cm  is  less  than  the  assumed  probe  radius  of  0.56  m,  the  simple  Lam  theory 
obviously  cannot  be  applied  at  this  higher  density.  The  difficulty  is  in  the  assumed  radius 
of  the  spherical  probe  corresponding  to  the  same  area  as  the  actual  cylindrical  probe.  At 
the  high  densities  where  Tq  is  quite  small,  it  would  be  more  logical  to  use  an  effective 
spherical  probe  with  a  radius  closer  to  the  radius  of  the  cylindrical  probe  (0*19  m).  Using  a 
radius  of  0.19  m  for  the  spherical  probe  at  a  density  of  10^/cm^#  Lam's  theory  yields  a 
potential  of  about  12  V,  which  is  considerably  larger  than  the  experimental  values.  This 
calculation  just  illustrates  the  rather  obvious  conclusion  that,  when  the  sheath  radius  is 
not  too  much  larger  than  the  probe  dimensions,  the  detailed  shape  of  the  probe  becomes 
more  important  and  one  should  use  a  more  correct  geometrical  model  of  the  probe,  if 
possible.  Even  at  the  lower  densities,  a  smaller  effective  radius  for  the  probe  would 
increase  the  calculated  potentials,  bringing  them  into  closer  agreement  with  the  measured 
values. 


Figure  6.  Probe  potential  vs.  electron  beam  current  for  different 
plasma  environments 
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Figure  7.  Probe  potential  v$.  ion  beam  current  for  different 
plasma  environments 

The  relative  potential  as  a  function  of  radial  distance  is  shown  in  Figure  8.  The 
potential  is  set  equal  to  zero  at  t/tq  =  1,  that  is,  when  r  is  at  the  virtual  cathode  located 
at  r  =  rQ*  The  potential  rises  until  the  probe  is  reached  at  r  =  R. 

Some  additional  curves  that  can  be  useful  for  planning  probe  experiments  are  given  in 
Figure  9.  One  curve  shows  the  mean  ionization  path  length  (A)  versus  the  density  (N)  of 
the  neutral  background  air  molecules  for  the  peak  ionization  cross  section  (assumed  to  be 
2  x  1Cf1&  cm2)-  The  significance  of  this  parameter  is  discussed  in  Section  3,4.2.  The 
other  set  of  curves  shows  the  radius  (tq)  of  the  sheath  region  as  a  function  of  the  bulk 
plasma  density  (n)  and  the  emission  current  of  electrons  (I)  for  a  plasma  electron 
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temperature  of  1000°  K.  For  ion  emission,  the  rQ  values  in  Figure  9  should  be  multiplied 
by  .  The  equation  shown  on  the  figure  for  <fu  .  is  the  probe  potential  at  which 

orbital  limiting  starts  to  occur.  It  is  obtained  from  Eq.  3-21  in  Section  3.5  by  solving  for 
♦,  letting  p  =  rQ,  and  identifying  1/2  mVQ2  with  kT.  For  values  of  probe 
potential  ^  >  4>q  ^  >  all  electrons  that  start  at  Tq  with  a  component  of  velocity  pointing 
toward  the  probe  will  be  captured  by  the  probe  and  thus  there  is  no  orbital  limiting. 
However,  for  l  ,  some  of  the  electrons  that  start  at  rQ  will  circle  about  the 

probe  without  striking  it,  which  means  orbital  limiting  does  occur.  This  figure  and 
equation  are  applicable  even  if  ionization  of  the  neutral  atoms  occurs  and  reduces  the 
probe  potential  below  the  Lam  value  for  no  ionization,  as  discussed  in  Section  3.4.  A 
similar  set  of  curves  that  is  only  applicable  without  ionization  is  shown  in  Figure  10.  In 
this  f  igure,  at  the  probe  potential  corresponding  to  a  given  current  and  Tq/R,  if  rQ/R  is  less 
than  the  impact  parameter  p/R,  orbital  limiting  will  not  occur.  For  example,  for  I  =  0.01  A 
and  (too  =  ^000  V,  r^/R  =  20.  If  kT  »  1  eV,  p/R  =  30,  so  orbital  limiting  does  not  occur.  On 
the  other  hand,  for  kT  =  10  eV,  p/R  =  10,  so  orbital  limiting  would  occur. 


Figure  8.  Potential  spatial  variation  away  from  probe 
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SHEATH  RAO I US,  rQ,  AND 
MEAN  IONIZATION  LENGTH.  X 


Ml  lltli/M  III  NM  IV  .  M  (mu  ') 

IU1  1  III1''  Ml1  1  Hi1 


Figure  9*  Sheath  radius  for  electron  beam  emission  and  ionization 
mean  free  distance 


no  air  ionization 
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3.4  COMPUTER  MODEL  AND  RESULTS 


In  the  following  discussion,  the  emitted,  and  therefore  the  attracted  particles  are 
taken  to  be  electrons  because  ionization  effects  will  normally  be  important  only  for 
electron  emission  and  the  main  purpose  of  this  code  is  to  study  ionization  effects, 

3*4.1  Model  Without  Ionization 

The  computer  model  without  ionization  is  very  similar  to  Lam's  model  in  the  depletion 
region.  The  only  differences  are  in  the  outer  boundary  condition  and  a  slight  difference  in 
the  equation  for  the  electron  velocity. 

In  this  model,  the  outer  boundary  of  the  problem  is  taken  to  be  the  edge  of  the 
depletion  region.  Instead  of  having  a  smooth  transition  from  the  depletion  region  to  the 
bulk,  as  in  Lam's  model,  the  present  model  assumes  a  sharp  transition  from  the  depletion 
region,  where  the  ions  are  completely  swept  out,  to  the  bulk  region  where  the  ion  density 
is  equal  to  the  bulk  density.  This  approximation,  which  is  the  same  one  used  for  the  series 
solutions  in  Section  3.2,  is  justified  for  relatively  high  electron  emission  currents,  where 
the  potentials  are  relatively  large  compared  to  the  ion  thermal  energy  and  the  sheath 
region  is  large  compared  to  the  transition  region.  Since  the  net  charge  outside  the 
depletion  region  boundary  is  assumed  to  be  identically  zero,  the  electric  field  at  the 
boundary  is  zero,  and  the  potential  $  is  defined  to  be  zero  at  that  point.  A  central 
element  of  this  model  is  that  the  outer  edge  of  the  sheath  is  given  by  the  condition  that 
the  inward  electron  thermal  current  just  equals  the  electron  emission  current  (Eq,  3-13). 

The  second  difference  between  the  computer  model  and  Lam’s  model  and  the  series 
solutions  is  that  the  radial  electron  velocity  av.  a  position  where  the  potential  is  $  is  given 
in  the  code  by 

v  =  /v*  ♦  2e<J>/m  .  (3-18) 

“  t  n 

where  e  is  the  absolute  value  of  the  electronic  charge.  This  difference  from  the  electron 
velocity  in  Eq.  3-4  is  significant  only  close  to  rQ  where  is  very  small.  If  Eq.  3-4  were 
used  instead  of  Eq.  3-18,  the  carrier  densities  (H/4Trr  ve,  Eq.  3-2)  close  to  r^  would  be 
considerably  greater  than  the  bulk  density  nQ,  and  the  average  radial  velocity  ve  would  be 
correspondingly  less  than  vt^.  The  use  of  Eq.  3-18  ensures  a  smooth  variation  of  nfi  across 
the  boundary  region.  No  analysis  has  been  performed  to  determine  how  much  this  change 
affects  the  overall  probe  potential.  However,  since  it  was  convenient  and  easy  to 
incorporate  into  the  code,  this  change  has  lieen  utilized. 

With  no  ionization,  it  is  a  straightforward  matter  in  the  code  to  integrate  Poisson's 
equation  inward,  starting  from  the  sheath  boundary  r  *  rg,  where  $  and  the  electric  field 
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both  equal  zero  in  this  model.  Using  a  finite-difference  mesh  grid  where  the  potentials, 
and  consequently  the  velocities  and  carrier  densities,  are  defined  on  the  mesh  points  and 
the  electric  fields  are  defined  at  the  half-mesh  points,  the  integration  can  proceed  from 
mesh  to  adjacent  mesh,  from  r  =  rQ  to  the  probe  radius,  r  =  R. 

Table  3  shows  a  comparison  of  calculated  probe  potentials  obtained  with  this  code 
(with  no  ionization  and  with  fully  depleted  ions  in  the  sheath  region)  with  results  from 
Lam.  The  agreement  is  quite  satisfactory  considering  the  differences  in  the  two  models. 


Table  3.  Comparison  of  Results  from  Computer  Model  with  Lam 


Plasma  Density 

-  103/cm3 

Electron  Thermal 

E  mission  Current 

Potential  (V)  for  1-m-radius  probe 

Energy  (eV) 

(A) 

L  am 

Computer  Model 
(No  Ions  or  ionizaton) 

10.0 

6  x  10~3 

28 

28.5 

10.0 

6  x  10'2 

852 

960 

0.1 

6  x  10'3 

184 

230 

With  the  computer  model,  it  is  also  possible  to  simulate  the  condition  of  electron 
emission  and  an  electron  return  current  with  a  fixed  background  of  positive  charge  equal 
to  the  bulk  ion  density.  This  calculation  would  correspond  to  the  early-time  ABORC 
calculations  before  the  ions  have  had  time  to  be  swept  out  of  the  sheath  region.  Table  4 
shows  a  comparison  of  the  probe  potentials  from  this  code  for  fixed  and  swept-out  ions  and 
with  the  corresponding  ABORC  calculations  with  fixed  ions.  The  comparison  of  the  fixed- 
ion  results  with  ABORC  is  quite  good  for  large  emission  currents,  especially  considering 
the  differences  in  the  geometry  and  calculat  ional  techniques.  The  reason  for  the 
differences  at  low  beam  currents  is  not  known  for  sure,  but  one  possibility  is  orbital 
limiting,  which  is  present  in  ABORC  but  not  in  the  steady-state  code.  From  Figure  10,  for 
1  =  6  x  10~3  A  and  <p  between  7  and  30  V,  the  curve  for  Tq/R  is  above  the  curve  of  p/R  for 
kT  =  10  eV.  Hence  orbital  limiting  is  occurring  for  this  case  and  the  potentials  should  be 
larger  than  predicted  by  the  steady  state  theory.  Therefore,  it  is  reasonable  that  the 
ABORC  code  should  give  a  larger  potential  for  this  case.  As  discussed  in  Section  3.5,  the 
orbital-limiting  effect  is  largest  when  the  electron  thermal  energy  is  relatively  large 
compared  to  the  probe  potential,  as  is  the  case  for  the  lower  emission  currents  in 
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Table  4.  For  more  realistic  thermal  energies  (~  0.1  eV),  this  effect  should  be  less 
important  even  for  probe  potentials  as  low  as  10  eV. 


Table  4*  Effect  of  Stationary  Ions 


(Plasma  Density  •  10 

W> 

Probe  Potential  (V)  ] 

— r— 

Electron 

Emission 

ABORC, 

Thermal 

Current 

Mjrrr 

1-m-radiu»  cylinder 

Energy  (eV) 

(A) 

■kQHI 

■BOB 

Fixed  Ions 

10 

OEM 

28.5 

mm 

30 

10 

960 

■9 

260 

10 

1 

15  kV 

2.7  kV 

0.1 

6  x  10'3 

230 

n 

— 

It  is  interesting  that  the  predicted  potentials  with  the  fixed  ions  are  considerably 
smaller  than  the  predicted  potentials  when  the  ions  are  assumed  to  be  completely  swept 
out  of  the  sheath  region.  The  physical  reason  for  these  lower  potentials  is  fairly  clear 
from  the  code  calculations.  When  the  ions  are  not  present,  the  electron  density  tends  to 
be  less  than  in  the  bulk,  in  spite  of  the  1/r  geometry  effect,  as  the  probe  potential 
accelerates  the  electrons  toward  the  probe.  On  the  other  hand,  with  the  ions  present,  the 
electron  density  tends  to  stay  roughly  equal  to  the  ion  density.  Therefore,  to  supply  the 
same  return  current,  the  high-density  electrons  in  the  fixed-ion  case  do  not  have  to  move 
as  rapidly  as  the  low-density  electrons  in  the  fully  depleted  case,  and  thus,  the  probe 
potential  can  be  smaller.  This  difference  in  potential  should  be  observable  in  a  time- 
dependent  calculation,  and  perhaps  in  an  experiment  with  a  low  emission  current  and  small 
probe  potential  such  that  ionization  is  not  a  major  effect.  When  the  electron  beam  is  first 
turned  on,  the  probe  potential  should  quickly  rise  to  a  plateau  corresponding  to  the  fixed- 
ion  case.  Then,  on  a  time  scale  corresponding  to  the  ion  sweepout  time,  the  probe  poten¬ 
tial  should  gradually  rise  toward  the  value  in  the  fully  depleted  case.  However,  this  time- 
dependent  effect  might  be  obscured  when  ionization  effects  are  large  and  have  a  major 
influence  on  the  probe  potential,  as  discussed  in  the  next  section. 


The  important  result  from  this  analysis  is  that  the  predicted  potentials  in  Table  3  with 
movable  ions  are  considerably  larger  than  the  measured  probe  potentials  with  electron 
emission  (Table  1).  Therefore,  some  important  physics  is  apparently  missing  from  the  non¬ 
ionization  model  used  thus  far* 

3.4*2  Model  with  Ionization 

When  an  electron  in  air  achieves  a  kinetic  energy  greater  than  the  ionization  energy 
of  the  air  molecules,  there  is  a  finite  probability  that  it  will  produce  additional  electron- 
ion  pairs.  If  this  ionization  occurs  in  a  region  where  an  electric  field  exists,  the  newly 
created  electrons  and  positive  ions  will  be  accelerated  in  opposite  directions,  and  thus, 
both  species  will  contribute  to  the  electric  current.  If  the  electric  fields  are  large  enough 
over  a  sufficient  distance,  the  new  electrons  will  reach  a  kinetic  energy  greater  than  the 
ionization  energy  and,  thus,  can  produce  additional  ionization.  This  multiplying  effect  is 
known  as  avalanching. 

To  simulate  ionization  effects  in  the  steady-state  computer  code,  the  ionization 
probability  curve  versus  electron  energy  was  approximated  as  a  linearly  increasing  curve 
starting  at  50  eV  and  going  to  200  eV.  Above  200  eV,  the  probability  curve  was  assumed  to 
be  fiat  for  all  higher  energies.  To  smooth  the  calculations  near  the  discontinuities  at  50 
and  200  eV,  a  smooth  parabolic  transition  was  used,  extending  10  eV  on  either  side  of  each 
discontinuity.  This  assumed  ionization  curve  has  a  somewhat  higher  threshold  energy  than 
the  experimental  curve  for  oxygen  in  Figure  14  (see  Section  4).  It  is  believed  that  the 
detailed  shape  of  this  ionization  curve  will  not  have  a  major  effect  on  the  computed 
results.  However,  this  shape  could  be  easily  modified  for  other  calculations  if  desired. 
The  peak  magnitude  of  the  probability  curve  is  an  input  parameter  in  the  code 
corresponding  to  different  densities  of  the  neutral  background  gas.  This  parameter  is 
given  in  terms  of  the  mean  distance  for  ionization  (  A)  by  an  electron  with  energy  greater 
than  200  eV  —  that  is,  with  the  maximum  ionization  probability. 

Ionization  can  come  from  three  different  sources  of  electrons:  (1)  emitted  electrons, 
(2)  return  electrons  which  originate  at  the  outer  boundary  of  the  sheath  region,  and  (3) 
electrons  created  by  the  ionization  which  are  then  accelerated  to  energies  above  the 
ionization  threshold.  In  the  present  code,  ionization  by  the  emitted  electrons  is  ignored 
because  their  density,  when  averaged  over  a  spherical  shell  around  the  probe,  is 
considerably  less  than  the  density  of  the  returning  electrons.  Also,  the  energy  of  the 
emitted  electrons  is  often  quite  high,  and  the  ionization  probability  actually  decreases  at 
high  energies,  contrary  to  what  was  assumed  above  for  the  code. 
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When  the  code  was  first  being  written,  it  was  intended  to  include  ionization  by  the 
newly  created  electrons  which  are  then  accelerated  to  the  ionization  threshold.  However, 
some  calculations  with  an  earlier,  simpler  version  of  the  code  indicated  that,  whenever 
ionization  is  significant,  the  potential  curve  becomes  relatively  flat  and  the  ionization 
electrons  seldom  reach  the  energy  for  the  ionization  threshold.  Therefore,  ionization  by 
these  electrons,  that  is,  avalanching,  has  been  omitted  from  the  code  for  simplicity, 
although  it  could  be  added  with  not  too  much  complexity  if  it  were  thought  to  be 
important*  Thus,  only  ionization  by  the  electrons  which  return  from  the  outer  boundary  of 
the  sheath  are  included  in  the  code. 

Ionization  is  incorporated  into  the  code  as  a  carrier  generation  rate,  dn/dt,  in  each 
mesh  region  proportional  to  the  density  of  electrons  passing  through  that  zone  and  as  a 
function  of  their  local  energy.  Both  the  electrons  and  the  ions  are  assumed  to  be  created 
with  zero  initial  velocities.  In  steady  state,  dn/dt  at  a  specific  mesh  station  j  will 
contribute  to  the  electron  current  and  density  at  all  mesh  stations  i  closer  to  the  probe 
than  j,  and  to  the  ion  current  and  density  at  all  mesh  stations  k  further  from  the  probe 
than  j. 

The  velocity  of  the  electrons  from  region  j  when  they  reach  region  i  is  taken  to  be 

ve(i)»  14.(1)  -  4>(  j  )  1  ,  (3-19) 

e 

where  <X0  and  M)  are  the  potentials  at  regions  i  and  j.  (In  the  sign  convention  used  in  the 
code,  the  potentials  are  positive  for  electron  emission.)  Since  physical  reasoning  indicates 
that  the  potential  curve  must  increase  monotonically  from  the  outer  edge  of  the  sheath  to 
the  probe,  the  argument  of  the  square  root  in  Eq*  3-19  is  always  positive. 

If  there  were  no  sweepout  velocity  due  to  the  motion  of  the  probe  relative  to  the 
ambient  plasma,  the  velocity  of  the  ions  from  region  j  when  they  reach  region  k  would  be 
given  by  an  equation  similar  to  Eq.  3-19  but  with  $0  replaced  by  <Kj),  <Ki)  replaced  by  <j>(k), 
and  me  replaced  by  the  ion  mass  m..  It  is  not  possible  to  simulate  rigorously  a  linear 
sweepout  velocity  in  a  spherically  symmetric  geometry.  However,  in  an  admittedly  rather 
crude  attempt  to  include  some  effect  of  sweepout,  the  ions  were  given  an  additional  radial 
outward  velocity  vs*  Thus,  the  ion  velocity  at  position  k  was  calculated  from  the  formula 

Vj(k)=  / vj  ♦  ^  t4>(i)  -  4>(k)l  •  (3-20) 

i 
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For  all  of  the  results  presented  later,  the  magnitude  of  mj  was  taken  to  be  2  x  10  me 
and  vs  was  taken  to  be  5000  m/s,  which  is  probably  somewhat  high  for  typical  probe 
velocities.  However,  again  it  is  felt  that  this  value  would  not  have  a  major  effect  on  the 
conclusions  from  this  study.  If  desired,  the  magnitude  of  v$  could  be  changed  for  future 
calculations. 

The  charge  density  p  in  Poisson's  equation  (Eq.  3-1)  is  the  algebraic  sum  at  each  mesh 
position  of  the  densities  from  the  inward-bound  electrons  that  originally  crossed  the  outer 
boundary  of  the  sheath  region  (as  in  the  zero-ionization  case)  and  the  electrons  and  ions 
created  by  the  ionization  process. 

The  steady-state  solution  to  this  problem  with  ionization  is  mathematically  very 
sensitive  at  high  emission  currents  and  relatively  small  mean  ionization 
distances  (  A  <  100  m).  The  reason  is  that  the  potential  curve  with  heavy  ionization 
becomes  very  flat  after  the  potential  reaches  the  ionization  threshold  around  50  eV. 
Therefore,  the  velocities  of  the  electrons  and  ions  that  were  created  by  the  ionization  due 
to  the  returning  electrons  (Eqs.  3-19  and  3-20)  become  quite  small  and  the  corresponding 
densities  become  very  large.  To  solve  this  problem  even  for  values  of  A  as  large  as  1000 
m,  it  was  necessary  to  use  a  calculat ional  approach  which  inverts  a  linearized  2N  x  2N 
matrix,  where  N  is  the  number  of  mesh  regions  used  in  the  calculations,  typically  about 
100.  The  factor  of  2  occurs  because  both  $  and  its  derivative  (E)  are  used  as  variables. 
Since  the  equations  are  very  nonlinear,  the  procedure  used  is  to  estimate  initial  potential 
and  electric  field  curves  over  all  of  the  meshes.  The  equations  are  then  linearized  for  per¬ 
turbations  in  <J>  and  E  from  these  initial  curves.  This  linearized  matrix  is  then  inverted  to 
obtain  the  changes  in  E  and 

In  principle,  one  would  like  to  use  these  new  curves  for  <p  and  E  around  which  to  again 
expand  the  equations  and,  thus,  repeat  the  process  until  the  calculated  changes  in  E  and  <f> 
are  less  than  some  selected  limit.  Unfortunately,  the  new  calculated  <f>  curve  is  often  not 
monotonically  increasing  across  the  sheath  region.  Consequently,  this  curve  cannot  be 
used  since  Eq.  3-19,  and  perhaps  Eq.  3-20,  would  give  imaginary  velocities.  Thus,  after  the 
linearized  changes  to  E  and  <j>  have  been  calculated,  some  other  algorithm  is  needed  to 
choose  a  new  estimated  <p  curve  that  is  monotonically  increasing.  There  is  an  infinite 
number  of  ways  this  could  be  done,  since  any  procedure  that  is  convergent  should 
theoretically  arrive  at  essentially  the  same  converged  result.  The  procedure  that  has  been 
used  is  to  adjust  the  old  values  of  E  by  some  fraction  FI  of  the  calculated  perturbations  in 
E  at  every  mesh  station,  provided  that  the  new  adjusted  E  is  not  negative.  If  any  new  E 
were  negative,  the  old  (positive)  E  at  that  mesh  is  reduced  by  some  fraction  F2  of  its 
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value.  These  estimated  E's  are  then  integrated  from  Tq  inward  to  obtain  the  corresponding 
potential  curve  about  which  the  equations  are  again  expanded.  There  is  no  guarantee  that 
this  procedure  will  always  be  stable,  especially  if  FI  and/or  F2  are  made  too  large  in  an 
attempt  to  speed  up  the  convergence.  In  fact,  at  large  currents  and  small  A's,  it  was  often 
difficult  to  obtain  satisfactorily  converged  solutions.  If  this  code  were  to  be  used  in  the 
future  for  those  more  difficult  situations,  this  adjustment  algorithm  should  be  reviewed 
and  perhaps  modified. 

All  the  calculations  with  this  code  thus  far  have  been  for  an  ambient  ionized  plasma 

i  3 

density  of  10  /cm  and  an  electron  temperature  of  1000°K  (kT  =  0.086  eV).  Figure  11 

shows  curves  of  the  calculated  potential  versus  distance  from  the  probe  for  an  emission 

current  of  0.1  A,  with  no  ionization  (  A  =  «)  and  A=  1000  and  100  m.  Note  in  particular 

the  flatness  of  the  potential  curve  for  A=  100  m  and  how  far  the  plateau  extends  from  the 

“16  2 

probe.  Assuming  a  peak  value  for  the  ionization  cross  section  of  2  x  10  cm  ,  a  value  of 
A  =  100  m  corresponds  to  a  density  of  neutral  background  atoms  of  5  x  lO^Vcm^.  A 
similar  calculation  for  A=  10  m  did  not  completely  converge,  but  the  plateau  potential  was 
approximately  50  eV.  Based  on  experience  with  the  other  A‘s,  it  is  expected  that 
essentially  this  potential  would  extend  in  to  the  probe.  Obviously,  ionization  distances 
even  as  large  as  100  m  have  a  very  pronounced  effect  on  the  magnitude  of  the  probe 
potential  to  infinity.  It  also  affects  the  shape  of  the  potential  curve  throughout  the  sheath 
region,  and  especially  close  to  the  probe.  Thus,  if  the  technique  for  determining  the  probe 
potential  experimentally  only  measures  the  differential  potential  between  the  probe  and  a 
point  in  space  a  few  meters  away,  this  measured  differential  potential  might  be  only  a 
small  fraction  of  the  probe  potential  to  infinity. 

It  is  fairly  clear  why  the  potential  approaches  a  constant  value  equal  to  the  ionization 
threshold  energy  at  relatively  small  values  of  A.  The  incoming  electrons  will  create 
electron-ion  pairs  more  or  less  uniformly  throughout  the  volume  near  the  probe  where  the 
local  potential  is  greater  than  the  threshold  energy.  The  electric  field  in  this  region  will 
accelerate  the  electrons  toward  the  probe  and  the  ions  away  from  the  probe.  The  electric 
field  that  results  from  this  separation  of  the  negative  and  positive  charges  will  decrease 
the  electric  field  that  had  existed  before  the  separation.  For  large  A  and  small  beam 
currents,  the  density  of  the  electron-ion  pairs  will  be  small  and  the  reduction  which  they 
produce  in  the  electric  fields,  and  therefore  the  probe  potential,  will  be  small. 
Conversely,  for  small  A  and  large  currents,  the  electric  field  produced  by  the  electron-ion 
pairs  can  almost  completely  cancel  the  previously  existing  field.  Of  course,  the  field 
cannot  become  negative  in  steady  state  because  the  electrons  must  always  be  attracted 


31 


toward  the  probe.  In  the  limit  of  very  high  ionization  density,  the  electric  field  must  be 
just  sufficient  to  remove  the  electrons  from  the  ionization  region  at  the  same  rate  that 
the  ions  are  removed,  primarily  by  their  assumed  sweepout  velocity,  Only  a  fraction  of 
an  electronvolt  potential  difference  is  required  to  produce  an  electron  velocity  equal  to 
the  assumed  value  of  v$  =  5000  m/s.  Hence,  in  that  case,  the  potential  curve  in  the 
ionization  region  will  be  very  flat  and  close  to  the  ionization  threshold  energy. 


Figure  11.  Radial  variation  of  potential  for  different  ionization  rates 

Figure  12  shows  similar  potential  curves  for  different  emission  currents.  The 
interesting  point  from  this  figure  is  that  the  probe  potential  to  infinity  actually  goes 
through  a  maximum  as  a  function  of  the  emission  current.  This  maximum  is  illustrated  in 
Figure  13,  where  the  potentials  at  the  probe  in  Figure  12  are  plotted  versus  the  emission 
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current.  At  small  values  of  current,  the  probe  potential  is  less  than  the  assumed 
ionization  threshold  (40  eV),  so  there  is  no  ionization  and  the  calculated  probe  potentials 
agree  with  the  Langmuir-Child  prediction.  Perhaps  just  as  significant,  the  differential 
potential  between  the  probe  and  a  point  a  few  meters  from  the  probe  is  much  larger  at  a 
current  of  0.01  A  than  it  is  at  currents  of  0.04  A  and  above. 


METERS  FROM  CENTER  OF  PROBE 


Figure  12.  Potential  profiles  for  different  emission  currents 


Figure  13.  Probe  potential  versus  beam  current 


The  reason  for  the  peak  in  the  l-V  curve  is  not  as  clear  as  the  reason  for  the  plateau 
in  the  potential  curve.  However,  a  plausibility  argument  can  be  developed  by  considering 
the  curves  in  Figure  12.  For  a  current  of  about  0.003  A,  the  probe  potential  is  about 
100  V,  which  is  well  above  the  ionization  threshold.  However,  the  radial  distance  over 
which  the  potential  exceeds  the  ionization  threshold  is  fairly  small  and  not  many 
ionization  events  occur  in  that  distance.  On  the  other  hand,  when  the  current  is  over  an 
order  of  magnitude  larger  (-0.04  A),  the  potential  reaches  the  ionization  threshold 
about  10  m  from  the  probe.  Hence,  there  is  a  long  distance  over  which  ionization  will 
occur  and  there  is  over  an  order  of  magnitude  more  returning  electrons  to  create  the 
ionization.  This  abundance  of  secondary  electrons  causes  the  potential  to  remain  closer  to 
the  ionization  threshold,  that  is,  less  than  the  potential  at  a  smaller  current.  It  is 
interesting  that  the  peak  in  the  l-V  curve  will  probably  be  less  pronounced  for  very  large 
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vehicles.  For  example,  from  Figure  12,  if  the  probe  radius  had  been,  say,  10  m,  essentially 
the  same  potential  curves  outside  that  radius  would  occur  for  the  same  currents.  Although 
other  curves  at  higher  currents  are  needed  to  prove  the  point  conclusively,  it  appears  that 
the  maximum  potential  for  a  10  m-radius  sphere  might  be  only  about  b0  V. 

Although  it  is  unclear  what  densities  of  neutral  atoms,  and  therefore  Vs,  should  be 
associated  with  the  data  points  in  Table  1,  the  calculated  probe  potential  in  Figures  11  and 
12  for  Vs  on  the  order  of  100  m  or  less  are  in  much  better  agreement  with  measured  probe 
potentials  for  electron  emission  than  are  the  predictions  such  as  those  in  Table  3,  which  do 
not  include  ionization. 

3*5  ORBITAL  LIMITING 

At  this  time,  a  brief  discussion  of  orbital  limiting  is  useful. 

Consider  a  particle  with  mass  and  charge  q,  which  has  a  velocity  vQ  in  a  region  where 
the  electrostatic  potential  is  zero.  If  this  particle  is  attracted  by  a  central  force  field 
toward  a  probe  with  radius  R,  it  will  strike  the  probe  if  its  initial  impact  parameter  is  less 
than 


where  is  the  probe  potential  relative  to  infinity  (Ref.  5).  Note  that  is  positive  if  the 
attracted  charge  q  is  negative.  If  its  impact  parameter  is  greater  than  p,  the  particle  will 

miss  the  probe  and  continue  on  its  pcith  about  the  probe. 

2 

If  we  identify  (1/2  m  vQ  )  with  the  thermal  energy  of  electrons  at  infinity  (kT  =  0.086 
eV  for  the  problems  in  Section  3.4)  and  assume  a  probe  potential  of,  say,  200  V,  and  R  =  1 
m,  the  critical  impact  parameter  p  is  48  m.  By  comparison,  the  outer  radius  of  the  sheath 
r0  for  an  emission  current  of  0.1  A  is  31.8  m  for  the  plasma  conditions  used  in 
Section  3.4.  Thus,  every  particle  at  rQ  that  has  a  component  of  its  velocity  directed 
toward  the  probe  has  an  impact  parameter  smaller  than  p  =  48  m  and  would  be  captured  by 
the  probe.  Consequently,  using  the  free-fall  approximation  for  the  electron  velocities 
(Eqs.  3-18  and  3-19)  appears  justified  for  such  a  potential. 

From  Figures  11  and  12,  the  probe  potentials  with  very  high  ionization  might  be  as 
low  as  40  eV.  If  we  use  this  potential  in  Eq.  3-21  instead  of  the  200  V  assumed  above,  the 
critical  impact  parameter  would  be  only  18.5  m.  This  would  indicate  that  the  free-fall 
approximation  would  be  inaccurate  for  emission  currents  greater  than  about  0.035  A  (r0  > 
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18.5  rn).  However,  when  there  is  a  large  amount  of  ionization,  a  more  realistic  criterion 
for  the  critical  impact  parameter  might  be  to  use  an  effective  probe  radius  R  in  Eq.  3-21 
equal  to  the  location  of  the  knee  in  the  potential  curve  --  about  15  to  20  m  for  I  =  0.1  A 
(Figure  11).  Obviously,  this  approach  would  produce  an  extremely  large  allowable  impact 
parameter.  The  rationale  for  this  criterion  is  that  the  incoming  electrons  start  to  produce 
electron-ion  pairs  when  they  reach  the  ionization  threshold,  that  is,  at  the  knee  in  the 
potential  curve.  These  ionization  electrons  are  created  with  very  small  initial  velocities 
and  all  of  them  are  eventually  captured  by  the  probe.  They  can  supply  the  difference  in 
the  return  current  to  the  probe  with  only  a  negligible  change  in  potential  if  all  of  the 
primary  electrons  strike  the  probe  in  one  case  or  if  some  of  them  miss  the  probe  and  orbit 
about  it  in  another  case.  Hence,  even  if  some  orbit  limiting  does  occur,  it  will  have  only  a 
minor  effect  on  the  probe  potential  when  there  is  high  ionization. 

A  more  rigorous  discussion  of  orbital  limiting  should  consider  a  Maxwellian 
distribution  of  the  initial  plasma  electrons,  as  was  done  in  Reference  5.  On  page  130  of 
Reference  5,  it  is  shown  that  the  current  collected  by  a  spherical  probe  with  no  ionization 
is  very  nearly  the  thermal  electron  current  that  crosses  the  sheath  boundary  at  r  =  rQ 
whenever 


exp 


q  'I'o/kT 


is  small  compared  to  unity.  This  quantity  is  the  exponential  e  *  in  Equation  45  of 
Reference  5.  For  =  200  V,  r^  =  31.8  m,  and  T  =  1000°  K,  the  magnitude  of  the  above 
exponential  is  about  0.1,  so  the  criterion  is  reasonably  well  satisfied.  Actually,  with  no 
ionization  and  this  current,  <t>oo  would  be  much  greater  than  200  V#  so  the  criterion  would 
be  satisfied  even  better.  With  ionization,  the  effective  radius  R  of  the  probe  should  be 
increased,  for  the  reasons  discussed  previously.  This  would  make  the  absolute  magnitude 
of  the  argument  of  the  exponential  larger,  even  for  ^  -  50  V,  so  the  exponential  would  be 
less  than  0.1. 

In  summary,  for  the  range  of  parameters  used  thus  far  in  the  computer  model,  and 
probably  for  the  major  range  of  interest,  orbital  limiting  appears  to  be  an  insignificant 
effect,  and  the  free-fall  approximation  for  the  attracted  electrons  is  well  justified. 
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4*  ABORC  CALCULATIONS 


This  section  describes  calculations  of  electron-beam  emissions  performed  self- 
consistently  in  2-1/2  dimensions  using  the  ABORC  SGEMP  computer  code  (Ref.  9).  "Self- 
consistent"  in  this  context  means  that  electromagnetic  fields  can  alter  charged  particle 
trajectories,  thus  affecting  spatial  currents  and  the  fields  themselves.  Advantages  of  the 
numerical  approach  described  here  are  the  minimization  of  environmental  assumptions, 
capability  to  meet  multidimensional  geometries,  and  ease  of  operation  of  a  previously 
existing  tool.  The  code  employs  a  full  description  of  negative  and  positive  particles  which 
can  move  self-consistently  in  electromagnetic  fields  obtained  from  finite-differenced 
Maxwell's  equations,  and  computes  the  time  evolution  of  fields  and  currents  around  the 
vehicle.  This  allows  for  effects  of  orbital  limiting,  space-charge  limiting,  and  the 
geomagnetic  field  automatically  in  the  simulations.  The  approach  has  many  similarities  to 
that  of  Rothwell,  executed  previously  in  one  dimension  (Ref.  10).  A  test  problem 
compatible  with  both  ABORC  and  RothwelPs  code  was  executed  at  AFGL's  direction  and 
gave  similar  results,  as  discussed  below. 

Disadvantages  of  using  particle-pushing  codes  such  as  ABORC  are  typically  the  large 
computer  core  and  execution  time  required  to  obtain  stability  and  statistical  accuracy  in 
certain  parameter  ranges  of  interests  (Ref.  16),  particularly  late-time  conditions. 
Therefore,  the  results  obtained  for  the  present  problem  with  the  computer  code  are 
limited  to  a  small  fraction  of  the  relevant  parameter  space— in  particular,  relatively  low 
plasma  densities  and  large  plasma  temperatures,  and  to  relatively  early  times  in  the 
complex  evolution  of  the  probe/beam/environment  system.  However,  usage  of  the  code 
elucidates  complex  behavior  caused  by  geometry,  space-charge-limiting,  and  magnetic 
field  effects,  which  would  be  much  more  difficult  to  address  with  simpler  analytic 
models.  Also,  relevant  parameters  for  analytic  models  can  be  determined  by  exercising 
the  particle  model  for  variations  in  the  parameters.  Results  can  be  employed  in 
constructing  other  models  capable  of  solving  the  entire  parameter  range  of  interest. 

The  remainder  of  this  section  describes  the  computer  code  and  modifications  added  to 
simulate  the  background  plasma  and  the  earth's  magnetic  field.  The  rudiments  of  the 
model  are  described  and  exercized  for  simple  cases.  Predictions  of  the  effects  of  geomag¬ 
netic  field  and  background  ionization  on  the  probe  response  are  presented  for  wide  ranges 
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of  beam  currents*  Finally,  a  discussion  of  results  in  terms  of  the  available  experimental 
data  is  given,  and  implications  for  future  measurements  of  interest  are  considered*  Some 
of  the  initial  modeling  and  results  leading  to  the  calculations  reported  here  are  discussed 
in  Reference  17. 

4.1  DESCRIPTION  OF  THE  STANDARD  ABORC  CODE 

ABORC  is  designed  to  solve  Maxwell's  equations  by  direct-finite-differencing  in  gene¬ 
ralized  coordinates  for  axisymmetric  geometries.  Spatial  current  densities  are  obtained 
from  finite  'particles*  of  charge  which  are  followed  through  the  spatial  mesh  of  zones. 
Each  particle  represents  many  negative  or  positive  charges  and  is  acted  on  by  the  local 
electric  and  magnetic  fields  during  each  time  step.  Emission  of  arbitrary  energy,  angular, 
spatial,  and  time  distributions  of  currents  can  be  specified.  The  calculational  volume  may 
contain  either  conductors  or  vacuum,  with  variable  conductivity,  permittivity,  and 
permeability  (a,  e,  p  ). 

Boundary  conditions  for  the  ABORC  code  require  the  specification  of  an  outer, 
perfectly  conducting  cylinder.  Free-space  solutions  can  be  obtained  by  moving  the  outer 
boundary  out  so  that  there  are  many  plasma  Debye  lengths  between  the  probe  and  the 
outer  walls  (Refs.  17-18).  Finite  conductivities  can  be  specified  representing  imperfect 
conductors,  and  dielectric  structures  and  high-permeability  regions  may  be  treated. 
Backscattering  of  electrons  can  be  specified  where  charge  is  re-emitted  from  surfaces 
upon  contact. 

The  code  was  typically  dimensioned  for  6,000  spatial  zones  which  may  vary  in  size. 
Up  to  100  conducting  regions  for  specification  of  bodies  of  revolution  can  be  employed. 
Randomizing  techniques  were  employed  for  electron  emission  distributions.  As  many  as 
8,000  particles  were  tracked  during  any  given  time  step  for  as  many  time  steps  as 
desired.  The  code  is  written  in  FORTRAN,  and  typical  computer  times  vary  from  3  to  60 
min  on  the  CDC  7600  (or  approximately  5  times  longer  on  the  CDC  6600). 

4.2  MODIFICATIONS  FOR  THE  PLASMA-PROBE  CALCULATION 

The  plasma-probe  calculations  required  additions  to  the  standard  code  for 
specification  of  the  ambient  ionized  plasma,  the  magnetic  field,  and  ionization  of  the 
neutral  particles  by  electrons.  The  modeling  for  these  additions  is  described  here,  along 
with  modifications  permitting  use  of  the  system  on  the  AFGL  CDC-6600  computers. 
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The  background  plasma  is  modeled  by  a  large  number  of  negative  macroparticles 
representing  electrons,  and  stationary  positive  macroparticles  representing  ions*  The 
positive  macroparticles  are  a  direct:  and  computationally  free  result  of  the  fact  that  the 
code  basically  creates  a  quasi-neutral  plasma  at  time  zero  and  then  moves  the  electrons 
while  the  ions  remain  fixed.  The  initial  energy  distribution  of  the  electrons  is  rectangular 
with  midpoint  temperature  value  X  (interpreted  loosely  as  the  plasma  temperature)  and 
with  a  spread  equal  to  T;  that  is,  the  distribution  extends  from  T/2  to  3  T/2  with  a 
constant  magnitude*  The  initial  angular  distribution  is  isotropic,  and  the  number  density, 
n,  is  constant,  ignoring  the  steady-state  sheath  that  actually  develops  around  a  spacecraft* 
This  is  a  good  approximation  as  long  as  the  probe  potential  with  the  beam  on  is  large 
compared  to  the  plasma  temperature,  and  therefore  large  compared  to  the  probe  potential 
before  the  beam  is  turned  on.  The  plasma  is  constantly  exiting  and  entering  at  the  outer 
boundaries  of  the  simulation  volume*  Assuming  that  the  boundaries  are  far  enough  away 
from  the  probe,  a  constant  flux  of  electrons  (equal  tonev/2,  where  v  is  the  average 
normal  velocity)  strikes  the  walls*  The  plasma  loss  is  compensated  by  a  constant  inward 
emission  of  plasma  particles  from  the  outer  boundary,  corresponding  to  the  thermal 
motion  only. 

The  simplifying  assumptions  in  the  plasma  description  are  thought  to  be  justified  for 
the  present  calculations*  Ions  move  approximately  1/10  of  a  probe  dimension  in  the  first 
10  ys  calculated  here,  so  a  zero  velocity  for  the  ions  introduces  very  little  error.  A 
similar  argument  applies  to  a  zero  assumed  rocket  velocity*  The  uniform  energy 
distribution  assumption  has  not  been  tested  for  accuracy  in  the  present  work,  but  it  is  not 
thought  to  introduce  any  significant  error* 

Effects  of  the  geomagnetic  field  are  modeled  with  the  field  along  the  axis  of  the 
cylindrical  coordinate  system  due  to  rotational  symmetry  requirements  in  the  code*  The 
appropriate  terms  were  included  in  the  Maxwelf s  equations: 
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where  each  term  which  references  E  ,  H7,  or  H  has  been  newly  added  for  the  effects  of 
the  geomagnetic  fieid.  The  particle  equations  of  motion  become: 


v*  *  v  +  (E  ♦  v  -  v.B  ) 
z  z  my  z  r0  0  r  ' 


v+  =  v  ♦  ^  (E  +  v„B  -  V  B0)  , 
r  r  my  f  0z  z  0  ' 


v*  =  vQ  +  ^  (Eq  +  v  B  -  v  B  )  , 
0  0  my  0  zr  r  z  ' 


/  =  y  ♦  a£t_  (ez  (v*  ♦  vz)  ♦  Er  (V;  +  vr)  *  Ee  (Vq  ♦  v0)j  , 


2mc 


where  the  superscript  ♦  denotes  values  at  the  end  of  the  time  step,  and  the  newly  added 
terms  in  the  equations  which  treat  the  geomagnetic  effects  are  those  referencing  E_,  B  . 

O  2- 

and  Br#  Also,  the  substitution  B  =  pH  has  been  employed,  where  pis  the  permeability  of 
free  space.  More  detail  on  the  development  and  solution  of  these  equations  is  given  in 
Reference  9. 

Simulations  begin  with  a  constant  value  for  the  geomagnetic  field,  Be^,  everywhere 
in  the  grid  at  time  zero.  As  particles  gradually  pick  up  a  net  azimuthal  velocity,  a 
magnetic  field  is  generated  in  the  opposite  direction  to  Be*  The  opposing  field  is,  in 
general,  spatially  dependent.  Hence,  the  particle  motion  and  fields  are  calculated  "self- 
consistently.* 

Ionization  of  neutral  species  is  modeled  using  the  cross  section  shown  in  Figure  14 

(from  Ref.  19)#  A  single  curve  for  atomic  oxygen  represents  ionization  by  electrons  of  all 

species  present  in  the  ionosphere  in  the  present  model.  The  accuracy  of  the  cross  section 

is  within  a  factor  of  2  for  the  other  species  over  most  of  the  energy  range,  which  is 

sufficient  for  the  present  model  accuracy.  In  particular,  the  peak  cross  section  for 

which  is  the  other  dominant  species  at  the  altitudes  of  interest  for  sounding  rockets,  is 
“1  2 

about  2.5  x  10  cm  and  it  also  occurs  around  100  eV.  In  the  code,  the  oxygen  cross 
section  curve  is  approximated  by  a  curve  fit  to  Figure  14.  The  code  emits  new  particles 
representing  ionization  products  when  the  calculated  charge  buildup  in  a  zone  due  to 
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ionization  reaches  some  fraction  of  the  ambient  electron  density.  Specifically,  a  particle 
is  emitted  in  a  cell  due  to  ionization  when  the  following  is  true. 

N  ||  ^  o  v  d(?  dt  >  fn  , 

where  f  is  an  input  parameter  (typically  set  equal  to  1),  dn/d<¥  is  the  energy  spectrum  of 
electrons  in  each  zone,  n  is  the  total  number,  v  is  their  velocity,  and  N  is  the  neutral 
number  density.  The  emitted  electron  is  assumed  to  have  1  eV  initial  energy  and  random 
isotropic  direction,  and  leaves  behind  a  stationary  ion. 


•  N/gov  d<f  CALCULATED  AT  EACH  SPATIAL  POSITION  FOR  BEAM  AND  PLASMA  ELECTRONS 

•  SECONDARY  ELECTRON  PARTICLES  EMITTED  FROM  EACH  CELL  WHEN  THE  CHARGE  DENSITY 

OUE  TO  IONIZATION  BECOMES  COMPARABLE  TO  INITIAL  BACKGROUND  ELECTRON  DENSITY 

•  SECONDARIES  EMITTED  WITH  1-eV  ENERGY 

•  COLLISIONAL  ENERGY  LOSS  AND  RANDOMIZING  OF  DIRECTION  OF  IONIZING  PARTICLES 

NEGLECTED 


Figure  14.  Summary  of  treatment  of  ionization  of  neutral  particles 
by  electrons  in  the  ABORC  particle  model  for  dynamic 
probe  behavior 
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43  COMPUTER  MODEL  FOR  THE  PROBE,  BEAM,  AND  PLASMA 

The  A  BO  K  C  model  for  the  probe,  electron  beam,  and  background  plasma  is  illustrated 
in  Figure  15.  Major  parameters  defining  problem  components  are  listed  on  the  right.  The 
symbols  are  used  in  the  remainder  of  this  section  to  define  problem  input  parameters.  The 
probe  has  been  considered  to  be  a  simple  cylinder  primarily  for  purposes  of  comparing 
initial  results  of  the  calculations  with  simpler  spherical  models.  A  more  realistic  rocket¬ 
shaped  geometry  could  be  treated  by  the  code,  and  may  behave  somewhat  differently  in 
magnitude,  but  will  probably  show  the  same  qualitative  behavior  as  a  function  of  major 
parameters.  It  is  assumed  that  electrons  which  strike  the  probe  are  not  re-emitted 
because  secondary  electrons  generally  have  low  energies  (around  2  eV),  so  they  would  be 
immediately  attracted  back  to  the  vehicle  and  thus  produce  negligible  effects.  However, 
a  problem  at  higher  plasma  density  or  lower  beam  energy  might  require  the  inclusion  of 
secondary  electrons*  Treating  the  vehicle  as  stationary  should  introduce  minimal  errors 
for  cases  where  the  ionization  of  background  neutrals  is  neglected.  The  time  scale  to 
move  an  object  dimension  is  much  longer  than  the  plasma  frequency.  A  stationary  ion 
model  should  be  sufficient  for  the  same  reason. 

Beam  parameters  define  the  energy,  angular,  spatial,  and  temporal  characteristics 
assumed  in  the  analyses.  The  values  shown  are  intended  to  represent  a  recent  experiment 
by  Cohen  (Ref.  2)  within  the  ability  of  the  computer  model.  In  actuality,  the  beams  are 
highly  concentrated  over  a  small  area,  are  emitted  at  some  angle  to  the  geomagnetic 
field,  and  are  confined  to  a  smaller  angular  spread.  In  the  code,  the  beam  is  assumed  to  be 
emitted  over  a  much  larger  area  to  ensure  that  several  finite-difference  zones  are 
overlapped  by  it.  The  rotational  symmetry  requirement  allows  a  pencil  beam  along  the 
coordinate  system  axis  parallel  to  the  geomagnetic  field  only.  For  emission  at  an  angle  to 
the  field,  the  beam  must  necessarily  fan  out  like  a  cone,  as  illustrated.  This  restriction 
introduces  questionable  geometry  differences  between  the  experiments  and  the  analytical 
models,  but  it  does  not  inhibit  magnetic  field  effects  on  the  beam  transport.  The  beam 
angular  spread  is  small  in  the  code,  as  it  is  in  the  experiments.  The  pulse  shape  is  a  step 
function  at  time  zero  which  accurately  models  the  time  interval  treated  by  the  code. 

The  pertinent  plasma  quantities  are  the  initial  electron-ion  number  densities  and 
temperatures  and  the  neutral  species  densities.  The  values  shown  are  for  the  ranges 
considered  in  the  present  calculations.  The  10-eV  electron  temperature  is  artificially  high 
(1  eV  is  the  physically  reasonable  upper  limit),  but  it  was  employed  in  most  cases  as  an 
economy  measure  during  the  development  stages  of  the  modeling.  The  positive  ions  have 
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been  taken  to  be  stationary  during  the  early  time  evolution  ( ~2Q  ms)  due  to  their  relatively 
large  mass.  The  computer  code  could  treat  them  as  mobile,  but  they  were  omitted  for 
reasons  of  economy. 


INPUT  PARAMETERS 


RE  *0302 3 


PROBE 

SHAPE 

SIZE 

MATERIAL 
VELOCITY,  v 


CYLINDER 

LENGTH  =  DIAMETER  =  2  m 

PERFECT  CONDUCTOR, 

ZERO  CHARGE  ALBEDO 

0 


BEAM 

ENERGY ,  <? 

ENERGY  SPREAD 
PITCH  ANGLE.  0 
PITCH  ANGLE  SPREAD 
CURRENT,  I  (amp) 
AREA 

RISE  TIME 
PULSE  WIDTH 


4. S  keV 
1  keV 
0  and  70° 

10° 

10~3,  10~2,  10"1,  10°,  10 
3 
0 


PLASMA  ENVIRONMENT 

AMBIENT  ELECTRON  o  .  .  c 

DENSITY,  n  (cm‘3)  10%  10%  10%  10D 

TEMPERATURE,  T  (eV)  10 

NEUTRAL  DENSITY,  N  UP  TO  1013/cm3 

ION/ELECTRON 
MASS  RATIO 


GEOMAGNETIC  FIELD, 

Be  (W/m?)  0  and  5.8  x  10 


Figure  15.  ABORC  code  model  for  a  rocket-borne  electron  beam  emission 
into  a  neutral  plasma 


Electric  and  magnetic  field  solutions  were  obtained  by  finite-differencing  Maxwell's 
equations  on  a  grid  of  28  axial  by  18  radial  cells.  The  minimum  size  was  0.33  m  near  the 
vehicle,  and  increased  smoothly  to  1.6  m  at  the  outer  boundary.  This  grid  permits 
approximately  two  cells  per  Debye  length  near  the  vehicle  and  is  economical  enough  to 
allow  solutions  of  the  field  equations  over  the  long  simulation  times.  The  mesh  size  is 
slightly  larger  than  the  Debye  length  in  the  outer  reaches  of  the  volume,  but  it  is  thought 
that  this  will  have  little  effect  far  from  the  object.  The  primary  limitation  to  finite- 
differenced  Maxwell's  equations  occurs  when  small  zones  are  required.  Then  time  steps 
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must  be  reduced  and  computer  times  may  become  excessive.  The  present  simulation  with 
a  0.74-m  Debye  length  and  non-space-charge-limited  beam  permits  a  reasonably 
economical  calculation. 

The  background  plasma  is  represented  initially  by  typically  2400  particles  represent¬ 
ing  the  free  electrons.  Considering  the  volume  of  each  particle  (particles  are  actually 
annuli  in  two  dimensions),  we  have  approximately  six  per  Debye  volume  near  the  vehicle. 
This  number  is  adequate  for  the  close-in  representation  of  the  plasma.  We  have  employed 
considerably  fewer  particles  per  unit  volume  at  larger  distances  from  the  object  for 
economic  reasons. 

The  initial  condition  of  the  satellite-plasma  system  was  artificial  in  the  simulations. 
In  an  actual  experiment,  an  equilibrium  would  occur  in  which  a  sheath  would  exist  about 
the  vehicle.  Our  simplified  initial  condition,  in  which  the  vehicle  instantaneously  appears 
in  the  middle  of  a  boundless,  homogeneous  medium  and  simultaneously  begins  emitting 
electrons,  results  in  system  adjustments  both  to  the  presence  of  the  conducting  object  and 
to  the  fields  caused  by  the  energetic  particles.  The  system  is  unaffected  by  this  simplifi¬ 
cation  after  a  few  plasma  periods,  but  the  time  evolution  may  be  suspect  at  very  early 
times. 

Time  scales  for  the  simulation  are  determined  by  the  Courant  condition  for  fields,  and 
by  the  minimum  transit  time  across  a  Debye  length  for  particles.  In  a  typical  case,  the 
step  sizes  employed  are  0.7  ns  for  fields  and  20  ns  for  particles.  Even  with  the  resulting 
30:1  time  step  ratio  of  particles  to  fields,  the  caiculational  cost  is  dominated  by  the 
particle  position  updating. 

The  following  summary  highlights  the  ABORC  code  particle  model  used  for  early-time 
dynamic  probe  behavior. 

•  Rotationally  symmetric  probe  and  beam  modeled  in  2-1/2  dimensions. 

•  Time-dependent  Maxwell*s  equations  solved  self-consistent ly,  including  effects 
of  space-charge  limiting. 

•  Conducting  outer  boundary  — 

♦b  =  °* 

J  =  nev/2,  where  v  is  the  inward  component  of  the  velocity  from  the 
initial  energy  distribution  emitted  isotropically. 


Lorentz  force,  including  geomagnetic  field  acting  on  beam  and  plasma  particles. 


Ionization  of  neutrals* 


•  Stationary  ions. 

•  E  lectron  albedo  =  0. 

4-4  CALCULATED  BEHAVIOR  OF  THE  PROBE  FOR  SIMPLE  CONDITIONS 

Figure  16  shows  the  time  evolution  of  the  probe  potential  for  the  simplified  conditions 
of  no  geomagnetic  field  and  no  ionization  of  neutrals,  for  several  assumed  background 
electron  densities.  Notice  that  the  background  reduces  the  potential  compared  to  the  case 
with  no  plasma  (the  dashed  line),  and  that  the  reduction  is  greater  for  greater  electron 
densities.  The  plasma  is  able  to  provide  replacement  current  through  its  thermal  velocity 
more  effectively  for  higher  densities  until,  at  the  highest  density  value,  the  probe 
potential  is  essentially  the  same  as  the  electron  temperature.  However,  as  pointed  out  in 
Section  4.2,  the  accuracy  of  the  calculation  becomes  marginal  when  the  probe  potential  is 
not  significantly  larger  than  the  plasma  temperature. 
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Figure  16-  ABORC  particle  calculations  for  probe  potential  time 
evolution  for  different  electron  densities 


The  case  with  no  background  plasma  shows  potential  increasing  linearly  with  time, 
since  no  charge  flows  back  to  the  probe  to  neutralize  it.  If  the  calculation  had  been 
extended  to  late  times,  the  potential  would  have  risen  to  the  beam  energy  (~4.5  kV),  and 
all  further  beam  electrons  leaving  the  probe  would  return. 
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A  direct  comparison  between  ABORC  and  Kothwell's  code  (Ref.  10)  was  made  using 
the  loVcm^  density  case  of  Figure  16  and  a  uniform  emission  from  the  surface  of  the 
ABORC  model  and  no  geomagnetic  field  to  be  consistent  with  the  Rothwell  model.  The 
steady-state  potentials  from  the  two  codes  were  within  2  or  3  V,  and  the  oscillation 
periods  were  within  20%,  which  is  excellent  agreement.  In  addition  to  the  different 
geometries  between  the  ABORC  and  Rothwell  models,  there  were  several  other 
differences  such  as  somewhat  different  boundary  conditions  and  the  Rothwell  model  is 
quasi-static  whereas  ABORC  is  fully  dynamic.  The  test  problem  calculations  are 
described  in  more  detail  in  Reference  17. 

4.5  EFFECTS  OF  THE  GEOMAGNETIC  FIELD 

When  the  earth's  magnetic  field  is  included  in  'he  calculat ions,  it  affects  both  the 
beam  and  plasma  electrons  (plasma  ion  motion  perturbations  due  to  the  B-field  are  negli¬ 
gible  because  of  large  mass  and  low  velocities).  The  effects  are  summarized  in 
Figure  17.  Beam  electrons  are  bent  into  Larmour  orbits  (radius  RL)  if  they  cross  magnetic 
field  lines.  They  can  reach  a  maximum  distance  perpendicular  to  the  field  of  2  R  L,  but  can 
spiral  to  great  distances  from  the  probe  along  the  flux  lines.  The  effect  is  dependent  on 
the  magnitude  of  the  beam  current  in  that  it  becomes  insignificant  at  large  beam  currents 
where  the  probe  potential  always  approaches  the  beam  energy.  Plasma  electrons  are  also 
constrained  to  move  in  tight  orbits  (equal  to  0.2  m  for  the  presently  employed  artificially 
high  temperature  of  10  eV).  Since  the  electrons  cannot  cross  field  lines,  return  current 
comes  from  a  tube  equal  approximately  to  the  probe  diameter.  Higher  probe  potentials 
result  from  the  smaller  area  of  the  plasma  supplying  return  currents. 

The  increased  probe  potential  due  to  the  geomagnetic  field  is  illustrated  dramatically 
in  Figure  18A  for  a  low-current  electron  beam  emission  parallel  to  the  field.  The  lower 
two  graphs  in  the  figure  are  the  total  electron  currents  returning  to  the  side  and  bottom 
(Figures  18B  and  18C)  of  the  probe  due  to  the  plasma.  Although  it  is  not  shown,  an  amount 
of  current  comparable  to  that  on  the  bottom  of  the  cylinder  returns  to  the  probe  on  the 
top  of  the  cylinder.  The  current  striking  the  bottom  is  only  moderately  affected  by 
inclusion  of  the  earth's  field,  B^,  because  electrons  move  parallel  to  it  in  striking  that  sur¬ 
face.  Current  returning  to  the  side  is  dramatically  reduced,  however,  because  the  low- 
energy  plasma  electrons  cannot  cross  the  field  lines.  The  effect  on  the  probe  potential  is 
to  cause  it  to  increase  substantially  over  the  case  without  Be  (upper  figure)  because  the 
plasma  cannot  provide  return  current  fast  enough  to  hold  it  down.  This  effect, 


demonstrated  here  numerically,  was  predicted  analytically  in  1967  by  Parker  (Kef.  6).  It  is 
interesting  that  the  sum  of  the  quasi-steady  current  to  the  side  of  the  cylinder  with  Bg  = 
0  (  *  4  mA),  the  relative  current  to  the  bottom  of  the  cylinder  (  «1  mA),  and  an  additional 
estimated  1  mA  to  the  top  of  the  cylinder  equals  the  emission  current  (6  mA),  as  it  should 
in  quasi-steady  state. 


EFFECTS  ON  BEAM 
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e 
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Figure  1 7.  Effects  of  earth's  magnetic  field 

Figure  19  shows  the  effect  of  the  geomagnetic  field  on  probe  potential  for  a  wide 
range  of  electron  beam  currents.  The  probe  potential  was  computed  with  and  without  the 
field,  Be#  for  a  low-density  plasma.  The  results  plotted  are  the  maximum  values  observed 
during  the  first  15  ms  of  the  simulation.  The  potential  which  would  have  been  achieved  at 
15  ms  in  the  absence  of  a  background  plasma  is  also  plotted.  This  latter  potential  would 
continue  to  increase  until  it  reached  the  beam  energy.  The  presence  of  Be  increases  the 
potential  when  there  is  a  plasma  by  about  an  order  of  magnitude  at  the  smaller  currents. 
Note  that  the  potential  finally  reaches  the  beam  energy  (~5  keV)  with  or  without  the 
magnetic  field  for  high  current  conditions  due  to  the  dominance  of  the  electric  field 
limiting  of  the  beam  particles. 

The  effect  of  the  geomagnetic  field  on  probe  potential  was  also  considered  for 
different  beam  pitch  angles.  The  beam  was  emitted  parallel  to  the  field  (0°  pitch)  and 
almost  perpendicular  to  it  (70°  pitch),  as  illustrated  in  Figure  20.  The  beam  was  restricted 
to  an  8-m  radial  distance  from  the  probe  in  the  latter  case,  but  it  spiraled  upward 
unrestricted.  The  final  probe  potential  was  found  to  be  the  same  in  either  case  for  a  low 
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beam  current  (~10  *  A)  with  a  5-keV  beam.  Effects  of  the  pitch  angle  may  be  slightly 
greater  at  high  beam  currents,  but  the  parameter  does  not  appear  to  be  significant  with 
regard  to  interpreting  probe  potentials  based  on  this  simulation. 
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Figure  18.  Probe  potential  and  plasma  return  currents  calculated  with 
and  without  the  geomagnetic  field 


NO  PLASMA  (t 

NO  PLASMA 

<t  *  15  ps ) 


B  -  0.58 
gaussy 


»  0 

€ 


♦beam 

—  —  ^  ^ 


5  keV 

n  *  103/cm3 
T  -  10  eV 
N  *  0 


□  COHEN 

90  eV  e'  8EAM 
at  150  km 


RE-03027 


10" 1  10u 
BEAM  CURRENT  (A) 


Figure  19,  Maximum  probe  potential  during  the  first  15  |is  calculated 

with  and  without  the  earth's  magnetic  field  using  the  ABORC 
code  particle  model 
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Hgure  20.  Illustration  of  effects  of  beam  pitch  angle  on  system  behavior 


f  ffec.ts  of  the  geomagnetic  field  on  the  early-time  plasma-probe  behavior  are 
summarized  below, 

•  B-field  hinders  plasma  motion  perpendicular  to  f  ield  lines. 

•  Probe  potential  increases  as  much  as  10  times  over  value  obtained  without  B- 
field  modeled. 

•  Effect  on  probe  potential  is  significant  at  all  but  highest  beam  current  levels. 

•  B-field  reduces  effective  collection  area  of  probe,  suggesting  a  potential 
dependence  on  probe  area  normal  to  field  lines. 

•  Probe  potential  is  not  dependent  on  beam  pitch  angle. 

The  field  hinders  plasma  current  returning  to  the  vehicle,  causing  vastly  greater 
electric  fields  than  would  otherwise  occur.  The  effects  of  the  geomagnetic  field  require 
modeling  for  all  but  the  highly  space-charge-limited  beam  currents  to  accurately  predict 
the  potential.  This  behavior  suggests  that  an  interesting  experiment  would  be  to  measure 
the  time-dependent  potential  of  a  rocket  with  the  major  axis  first  parallel  and  then 
perpendicular  to  the  field  lines.  The  small  area  normal  to  the  geomagnetic  field  would 
presumably  result  in  a  potential  that  is  much  greater  than  in  the  orientation  with  the  large 
area  normal  to  the  field  for  the  same  beam  emission  current.  The  calculations  showed 
that  the  potential  was  not  strongly  dependent  on  the  beam  pitch  angle  to  the  magnetic 
field. 

4.6  EFFECTS  OF  IONIZATION  OF  NEUTRALS  ON  EARLY-TIME  BEHAVIOR 

Inclusion  of  ionization  of  the  neutral  particles  in  the  background  environment  allows 
beam  and  plasma  electrons  to  create  additional  electrons  and  positive  ions  through 
collisions.  Effects  increase  with  increasing  neutral  density  in  this  regime,  and  with 
decreasing  beam  energy.  The  degree  of  ionization  will  be  dependent  on  the  magnitude  of 
the  beam  current  in  that  the  beam  current  affects  the  probe  potential  and  the  energy  of 
the  returning  electrons,  and  thus  the  average  ionization  cross  section  of  the  electrons. 
The  effect  of  ionization  in  general  is  to  reduce  the  probe  potential  by  increasing  the 
supply  of  electrons  available  to  allow  replacement  of  the  beam  charge  at  lower  field 
levels. 

Background  ionization  effects  are  shown  in  Figure  21,  where  the  probe  potential  to 
infinity  is  compared  with  and  without  a  high-density  neutral  species  present.  Notice  that 


50 


the  potential  is  approaching  the  beam  energy  (<?=  5  keV)  in  the  case  without  the 
background  species,  and  that  it  reaches  only  H  kV  with  ionization*  The  time  at  which 
ionization  effects  become  observable  on  the  graph  M  ps)  corresponds  well  with  the 
"dissipation  time,"  =  1/N w,  where  a  is  the  cross  section  corresponding  to  elect rons  of 
velocity  v  which  ionize  the  neutral  background  species  of  density  N.  is  the  time 

required  to  produce  ionization  charge  equal  to  the  initial  ambient  value,  n.  The  values  for 
Tdis  can  ^  as  *ow  as  MS  *or  Present  conditions,  if  the  ionizing  particles  are  in  the 
100-eV  region*  This  energy  is  no  doubt  prominent  in  the  plasma  distribution  for  the 
exhibited  probe  potential.  Incoming  electrons  are  accelerated  into  this  energy  range  and 
efficiently  reproduce  themselves  through  ionization* 


Figure  21*  Probe  potential  calculated  with  and  without  ionization 
of  a  neutral  background  species 


Unfortunately,  the  calculation  had  to  be  stopped  shortly  after  the  peak  potential  was 
achieved  due  to  computer  costs*  However,  notice  in  Figure  21  that  the  potential  is  still 
falling  at  9  ps.  It  might  have  converged  to  a  steady-state  value  much  closer  to  the  value 
typically  observed,  100  V,  if  it  had  been  carried  to  later  times. 
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The  importance  of  ionization  on  the  early-time  probe  behavior  was  quantified  over  a 
wide  range  of  beam  currents  by  computing  the  potential  with  and  without  ionization  of  the 
neutral  particles  included.  Both  sets  of  calculations  included  the  geomagnetic  field 
treatment,  and  the  ionization  case  was  done  for  the  estimated  maximum  neutral-particle 
density  of  10  neutrals/cm  .  Only  moderate  (<  factor  of  3)  decrease  in  the  peak  probe 
potential  was  computed  due  to  ionization  for  the  conditions  shown  in  Figure  22  during  the 
15  us  of  the  simulations.  The  minimum  mean  free  path  for  ionizing  collisions  in  the 
calculations  was  ~6  m,  and  so  the  lack  of  effect  for  early  times  is  not  unexpected  since 
typical  electrons  travel  at  most  a  couple  of  mean  ionization  lengths  in  the  simulation 
time  ( *9  ms)  ,  so  not  many  ionization  events  occur.  Late-time  (~  ms)  effects  can  still  be 
significant,  however  (see  Section  3). 

The  effects  of  the  ionization  of  neutrals  on  early-time  probe  behavior  are  summarized 

3  3 

below  for  a  low-density  ambient  environment  (n  =  10  /cm  ). 

•  Ionization  by  primary  electron-beam  particles  is  negligible. 

•  Acceleration  of  plasma  electrons  to  energy  >15  eV  causes  beam-current- 

magnitude-dependent  ionization  of  neutrals. 

•  For  a  5-keV  beam  and  low-ambient-density  electron  environment,  maximum 

effect  is  seen  at  beam  current  magnitude  of  0.06  A. 

•  Peak  potential  is  reduced  to  as  little  as  40%  of  value  computed  without 
ionization  for  a  neutral  density  corresponding  to  the  highest  value  in  the 
altitude  range  of  interest. 

•  Effect  on  probe  potential  diminishes  as  beam  current  is  increased  beyond  0.06  A 

as  the  probe  potential  approaches  the  beam  energy.  For  a  larger  beam  energy, 

the  largest  effect  of  the  ionization  would  probably  occur  at  a  larger  beam 
current. 

•  Effects  may  be  larger  at  higher  ambient  electron  densities. 

13  3 

In  Figure  22,  the  curve  for  N  =  10  /cm  is  less  steep  than  the  curve  for  N  =  0,  until 
the  potentials  approach  the  beam  energy.  Thus  the  effect  of  ionization  apparently 
increases  with  increasing  beam  current.  Based  on  the  steady-state  results  in  Section  3, 
this  dependence  probably  occurs  because  the  region  over  which  the  potential  in  space 
exceeds  the  ionization  threshold  is  larger  at  larger  currents.  Consequently,  there  is  more 
volume  in  which  ionization  can  occur,  and  thus  ionization  should  have  a  larger  effect  on 
the  probe  potentials. 
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Figure  22.  Early-time  peak  probe  potential  calculated  with  the  ABORC 
particle  model  with  and  without  secondary-electron  pro¬ 
duction  due  to  ionization  of  neutral  species.  The  value 
plotted  is  the  maximum  value  observed  during  the  first 
15  |is  of  beam  emission. 


4.7  SUGGESTED  EXPERIMENTAL  MEASUREMENTS  BASED  ON  CALCULATED 

EARLY-TIME  PLASMA/PROBE  BEHAVIOR 

Results  of  the  ABORC  code  calculations  of  the  early-time  behavior  of  the  probe  are 
not  directly  comparable  to  available  experimental  data  because  the  experimental  data  are 
for  essentially  the  steady-state  late-time  regime.  However,  the  differences  between  the 
predicted  early-time  and  measured  late-time  results  suggests  some  additional 
experimental  measurements. 

The  probe  potential  was  tending  toward  measured  low  values  when  the  calculations 
were  stopped,  after  having  shown  potentials  generally  much  greater  than  the  values 
measured  experimentally  on  a  1-ms  time  scale  (Figure  21).  For  example,  a  beam  current 
of  0.8  A  resulted  in  a  3-kV  peak  probe  potential  with  ABORC  compared  to  experimentally 
measured  values  of  3  to  30  V  (see  PRECEDE  results  in  Table  1).  The  slight  differences  in 
beam  energy  between  experiment  and  calculation  would  not  account  for  any  major 
discrepancy  in  the  results.  An  interesting  measurement  then  would  be  to  obtain  the  peak 
probe  potential  as  well  as  its  ’steady-state*  value.  Then  it  would  be  significant  to  see  if 
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the  peak  value  is  sensitive  to  the  rocket  area  normal  to  the  geomagnetic  field,  and  if  the 
steady-state  value  is  not.  A  rationale  for  the  measurements  is  that  the  early-time 
potential  is  strongly  dependent  on  geomagnetic  field  effects  which  may  be  dependent  on 
the  rocket  orientation,  whereas  the  steady-state  value  is  driven  largely  by  the  ionization 
cross  section. 


5.  SUMMARY  OF  RESULTS 


( 


j 

i 
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The  following  are  the  main  conclusions  from  this  study. 

1.  For  small  magnitudes  of  the  beam  current,  the  steady  state  potentials  of  an  ion-  or 
electron-emitting  probe  in  a  partially  ionized  plasma  can  be  predicted  fairly  well 
using  one-dimensional  spherical  models  without  ionization  of  the  netural  background 
particles* 

2.  For  large  ambient  plasma  densities,  where  the  Debye  length  is  small  compared  to  the 
vehicle  dimensions,  care  must  be  used  in  choosing  the  effective  spherical  radius  of  the 
probe  for  making  the  estimates. 

3.  For  electron-emitting  probes,  when  the  probe  potential  exceeds  the  threshold  for 
ionization  of  the  neutral  background  gas  and  the  mean  ionization  distance  by 
electrons  is  less  than  about  200  m,  it  is  important  to  include  in  the  predictions  the 
effect  of  ionization  of  neutral  particles  by  the  electrons  which  return  to  the  probe  to 
replace  the  beam  current. 

4.  The  l-V  curve  for  an  electron-emitting  probe  when  ionization  is  important  apparently 
goes  through  a  maximum.  At  very  large  electron-beam  currents,  the  probe  potential 
asymptotically  approaches  the  ionization  threshold  energy. 

5.  When  ionization  of  the  neutral  background  gas  is  important,  the  spatial  variation  of 
the  potential  curve  in  the  region  around  the  probe  will  be  quite  flat  for  considerable 
distances  away  from  the  probe.  Thus,  a  measurement  of  the  differential  potential 
from  the  probe  to  some  nearby  point  in  space  could  be  considerably  less  than  the 
potential  of  the  probe  to  the  ambient  plasma. 

6.  At  early  times  after  the  emission  beam  has  been  turned  on,  the  probe  can  rise  to  a 
potential  that  is  considerably  larger  than  the  late-time  potential  if  ionization  of  the 
neutral  gas  is  significant. 

7.  The  presence  of  the  earth*s  magnetic  field  (Be)  increases  the  probe  potential  relative 
to  what  it  would  be  in  the  absence  of  Be#  The  magnitude  of  this  change  is  dependent 
on  the  beam  energy,  the  beam  current,  the  densities  of  the  ambient  plasma  and  the 
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neutral  particles,  and  probably  the  orientation  of  the  probe  geometry  relative  to  the 
direction  of  the  magnetic  field.  However,  the  orientation  angle  of  the  emission  beam 
relative  to  the  direction  of  the  magnetic  field  has  little  effect  on  the  probe  potential. 


i 
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APPENDIX  A 

TWO-FLUID  PLASMA  CODE 


A.1  INTRODUCTION 

This  appendix  describes  an  alternative  approach  that  might  be  used  in  future 
calculations  of  plasma  probe  problems.  The  idea  for  this  approach  arose  late  in  the 
present  program  when  the  difficulties  of  simulating  late-time  effects  with  particle¬ 
pushing  codes  became  apparent.  A  small  amount  of  work  and  code  development  were 
performed  on  this  approach  during  this  program,  and  a  few  preliminary  results  were 
obtained.  However,  the  program  came  to  an  end  before  the  code  could  be  fully  developed 
and  checked  out  for  the  conditions  appropriate  to  realistic  probe  environments.  Hence, 
this  appendix  just  summarizes  the  major  ideas  for  the  approach  and  the  status  of  the  code. 

This  two-fluid  plasma  code  is  a  modification  of  a  previously  existing  code,  PR  ECHO, 
to  consider  the  plasma  electrons  and  positive  ions  as  two  independent,  charged, 
compressible  fluids.  The  original  PRECHG  code  was  basically  a  Poisson's  solver  which 
calculates  the  potential  on  an  r-z  mesh  grid  in  a  cylindrical ly  symmetric  geometry  for  an 
arbitrary  distribution  of  charge  in  the  radial  (r)  and  axial  (z)  directions.  That  code  starts 
by  moving  charge  from  one  point  in  the  geometry  to  another  along  r  or  z  mesh  lines, 
leaving  behind  a  compensating  amount  of  charge  of  opposite  sign  on  the  original  sites.  To 
conserve  electric  flux,  an  initial  electric  field  is  specified  along  these  mesh  lines, 
corresponding  to  the  amount  of  charge  that  was  transported  along  them.  Although  this 
initial  distribution  of  charge  conserves  electric  flux,  it  normally  would  not  satisfy 
Poisson's  equation  because  V  x  E  would  not  be  zero  everywhere  in  the  volume.  The  code 
then  iterates  a  series  of  equations  alternately  in  the  r  and  z  mesh  rows,  which  tend  to 
reduce  V  x  E  to  zero  everywhere  in  the  volume  of  interest.  The  iteration  is  terminated 
when  the  maximum  absolute  value  of  V  x  E  at  any  mesh  location  is  below  some 
prespecified  magnitude. 

When  the  code  is  modified  for  a  two-species  plasma,  the  net  current  flow  (ions  minus 
electrons)  along  each  mesh  line  at  the  beginning  of  a  time  step  is  used  to  calculate  a 
change  in  the  electric  field  along  those  mesh  lines.  Again,  these  electric  fields  normally 
do  not  correspond  to  V  x  E  =  0  everywhere  in  the  volume,  so  the  code  iterates  the  fields 
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after  each  time  step  to  bring  V  x  E  back  within  acceptable  limits.  The  change  in  electron 
and  ion  densities  at  each  mesh  point  is  then  calculated  from  the  continuity  equation. 
Finally,  these  carrier  densities  and  electric  fields  are  used  to  calculate  new  current  densi¬ 
ties  for  the  start  of  the  next  time  step,  and  the  process  is  repeated. 

Since  the  above  procedure  solves  only  Poisson's  equation  and  not  Maxwell's,  the 
results  are  necessarily  quasi-static  but  still  time-dependent.  Fortunately,  for  many 
situations  the  frequencies  of  interest  are  low  enough  that  a  quasi-static  solution  is  quite 
adequate. 

The  main  uncertainty  in  the  problem,  as  discussed  later,  is  the  proper  physics  to  put 
into  the  code  to  generate  the  local  current  densities  (as  functions  of  the  densities  and 
electric  fields  or  potentials)  for  the  carrier  species  that  is  attracted  to  the  probe  and  the 
one  that  is  repelled  from  the  probe. 

The  first  assumption  that  was  used  corresponded  to  both  species  being  in  quasi- 
equilibrium,  corresponding  to  the  local  potential  at  every  position  in  space.  This 
assumption  is  appropriate  for  relatively  high  densities  of  the  charged  or  neutral  particles, 
such  that  mobility-limited  flow  is  applicable.  The  resulting  equations  are  identical  in  form 
to  the  equations  for  electrons  and  holes  in  semiconductors.  The  present  authors  have  had 
extensive  experience  in  programming  and  solving  these  equations  for  semiconductors, 
which  was  readily  adaptable  to  the  plasma  problem.  This  model  gave  very  encouraging 
comparisons  to  the  analysis  of  Reference  10  for  a  non-emitting  probe,  where  the  resulting 
probe  potentials  were  not  very  large,  even  though  the  analytic  model  in  Reference  10 
employed  the  orbit-limiting  effect  for  the  attracted  (ion)  species.  However,  it  was 
recognized  that,  for  lower  plasma  densities  and  high  emission  currents  (large  probe 
potentials),  the  quasi-equilibrium  approximation  for  the  attracted  species  could  be 
significantly  in  error.  The  "free-fall"  approximation  for  the  attracted  species  is  more 
suitable  for  those  conditions.  Under  this  assumption,  the  velocity  of  the  attracted  species 
in  steady  state  is  determined  uniquely  by  the  local  electrostatic  potential  relative  to  the 
background  plasma.  Unfortunately,  the  present  program  came  to  an  end  before  the 
physics  corresponding  to  this  free-fall  approximation  could  be  fully  incorporated  and 
checked  out  in  the  code.  However,  preliminary  results  indicate  that  there  will  apparently 
be  no  fundamental  difficulties  in  solving  the  problem  using  this  approximation. 

A *2  MESH  GEOMETRY 

The  mesh  geometry  for  this  code  is  illustrated  in  Figure  23.  A  body  (conducting  or 
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dielectric)  of  arbitrary  r-z  dimensions  is  enclosed  inside  an  outer  conducting  cylinder. 
This  type  of  system  is  often  called  2-1/2  dimensions  because  bodies  of  finite  size,  similar 
in  shape  to  real  space  vehicles,  can  be  simulated.  If  the  outer  cylinder  is  located  more 
than  a  couple  of  sheath  thicknesses  from  the  inner  body,  it  will  have  no  significant  effect 
on  the  solution  to  the  problem. 


z 


Figure  23.  Mesh  geometry  for  code  (see  discussion  in  Section  A«2) 


Carrier  densities  are  defined  at  the  intersections  of  the  mesh  lines  and,  therefore,  on 
the  boundary  of  the  outer  cylinder  and  the  body.  These  densities  are  essentially  assumed 
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uniform  over  the  rectangular  regions  around  each  mesh  intersection  illustrated  by  the 
dashed  lines  in  Figure  23.  Electric  fields  and  current  densities  are  defined  along  the  mesh 
lines  but  at  locations  midway  between  the  mesh  intersections  —  that  is,  midway  between 
the  carrier  densities.  In  this  system,  the  continuity  equation  expresses  the  rate  of  change 
of  the  average  density  inside  one  of  the  dashed  rectangles  in  the  figure  due  to  the  net 
current  that  crosses  its  four  boundaries.  Similarly,  Poisson's  equation  requires  that  the  net 
electric  flux  leaving  one  of  these  rectangles  must  correspond  to  the  net  charge  (ions  and 
electrons)  inside  the  rectangle.  The  requirement  that  V  x  E  =0  is  rigorously  equivalent  to 
requiring  that  the  line  integral  of  the  electric  field  around  any  enclosed  area  defined  by 
mesh  lines  (solid  lines  in  Figure  23)  should  be  zero. 


A3  PLASMA  EQUATIONS 

The  basic  equations  for  this  code  involve: 

1.  Particle  current  densities,  Jpi,  for  each  species  between  adjacent  mesh 

e,  i 

intersections  (e  =  electrons,  i  =  ions). 

2.  Continuity  equations  for  each  species  at  each  mesh  intersection. 


n  .  a  V 

e,  i 


e,  i 


(A-1) 


where  n  is  the  average  carrier  density  in  a  mesh  region. 


Time  rate  of  change  of  electric  field  between  adjacent  mesh  intersections  due 
to  current  densities,  before  adjusting  V  x  E  toward  zero  (essentially  a  one¬ 
dimensional  Poisson's  equation), 


E  =  (‘Vi  - 


( A-2) 


where  £q  is  the  permittivity  of  free  space  and  qj  and  qe  are  the  charges  on  the 
ions  and  electrons,  respectively  (qj  =  ♦,  qe  =  -).  When  q  is  used  without  a 
subscript,  it  denotes  the  absolute  value  of  the  electronic  charge. 

4.  V  x  E  iteration  to  reduce  the  curls  to  essentially  zero. 


AJ.1  Current  Equation  for  Repelled  Species 

In  the  literature  on  charged-particle  probes,  it  is  generally  assumed  that  the  plasma 
species  that  is  repelled  from  the  probe  is  essentially  in  thermal  equilibrium  with  a 
Maxwellian  energy  distribution  at  every  point  in  space  (Ref.  5).  If  the  species  is  actually 
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in  equilibrium  (no  currents  and/or  rates  of  change  of  density),  the  particle  density  at  any 
place  in  space  is  given  by  the  equation 


n(x)  -  nQ  exp[+q4<x)/kT]  ,  (A-3) 

where  n^  is  the  particle  density  for  that  species  in  the  plasma  far  from  the  probe  where  <f> 
is  zero,  <J<x)  is  the  electrostatic  potential  at  position  x  (relative  to  zero  at  infinity)  = 
-  ^  E(r)dr,  T  is  the  Maxwellian  temperature  which  characterizes  the  particular  species, 
and  k  is  Boltzman's  constant. 

If  there  are  any  currents  of  the  repelled  species  ir,  the  plasma,  either  steady-state  or 
transient,  Eq.  A-3  is  not  rigorously  satisfied.  However,  it  is  fairly  close  to  correct  and  can 
be  used  to  develop  an  equation  for  an  effective  particle  current  density  (])  in  terms  of  the 
gradients  of  the  particle  density  and  the  potential. 

Consider  a  very  simple  system  that  has  only  two  equal-size  mesh  zones,  as  illustrated 
in  Figure  24,  with  densities  2  defined  at  points  1  and  2  along  with  the  corresponding 
electrostatic  potentials  ^  2.  When  the  particles  are  in  nearly  thermal  equilibrium  at  each 
point,  the  velocities  of  the  particles  with  densities  n^  and  n2  are  distributed  isotropically 
with  Maxwellian  distributions.  Therefore,  half  of  the  particles  at  n-j  will  be  moving  toward 
n2  with  an  average  velocity  in  the  positive  x  direction  given  by 


f  (  m  >i1/2  ,  2  „  / 2kJ 

v  =  I  v  (— —  J  exp(-mv  /kT)  dv  =  / - 

T  *  x  K2  rrkT J  x  x  n  m  ' 

x  0 


(A-4) 


where  m  is  the  mass  of  the  particles  (Ref.  5).  Similarly,  half  of  the  particles  at  n2  will  be 
moving  toward  n^  with  the  same  average  velocity  vy^  .  However,  due  to  the  potentials  ^ 
and  <J>2,  the  Maxwellian  distributions  are  skewed  in  the  direction  of  the  electric  field  (x 
direction  in  the  present  problem).  The  differences  in  the  potentials  at  points  1  and  2  from 
the  value  at  the  midpoint  between  1  and  2  are  t E  Ax/2.  Therefore,  if  the  potential  in  the 
exponential  in  Eq.  A-4  is  adjusted  for  the  electrostatic  potentials  at  points  1  and  2,  the 
average  velocity  of  the  particles  in  one  region  toward  the  other  region  is 

v  =  vr  exp(qE  Ax/kT) 

1  x 

and 

v  *  vT  exp(-qE  Ax/kt)  .  (A-5) 

*2  '  x 

The  rates  of  change  of  the  densities  n^  and  n2  in  this  two-zone  system  are  then 


(A-6) 
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Expanding  n^  and  n2  about  their  average  value  lna ve  =  1/2(0^  ♦  n2)]  and  the  exponentials  in 
Eq.  A-5  up  to  first-order  terms,  and  dropping  the  higher-order  products  of  An  =  n^  -  n2 
and  qE  Ax/kT,  Eq.  A-6  becomes 
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( A-8) 


Comparing  Eq.  A-8  to  the  continuity  equation  (n  =  -V*)),  one  can  define  an  effective 
current  density  crossing  the  boundary  between  mesh  points  1  and  2  as 

Jeff  =_r  '  n2  4  nave  ^  *  (A_9) 


There  is  a  conceptual  difficulty  with  this  effective  current  density  if  one  lets  the  mesh 
width  Ax  approach  zero.  In  that  case,  ]  also  approaches  zero  (for  a  smoothly  varying 
density).  However,  the  important  point  to  remember  is  that  le^  is  used  in  the  code  only 
in  a  form  equivalent  to  Eq.  A-8,  which  is  independent  of  the  mesh  width,  as  it  should  be, 
for  a  smoothly  varying  density  and  an  approximately  constant  E  (that  is,  Ax  «  a  Debye 
length). 
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Figure  24*  Two-mesh  zone  problem 


It  will  be  noted  that  Eq.  A-9  has  the  same  form  as  the  current  equations  for  electrons 
and  holes  in  semiconductor  theory  if  vT  Ax/2(kT/q)  is  identified  as  an  effective  mobility. 
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This  similarity  in  the  equations  allows  use  of  many  of  the  analytical  and  computational 
techniques  that  have  been  developed  to  study  semiconductor  problems  (Ref.  21).  However, 
it  is  emphasized  that  this  similarity  is  purely  format,  and  it  is  not  implied  that  the  repelled 
plasma  has  an  actual  mobility  nor  that  this  formulation  is  valid  only  in  a  very  high-density 
plasma.  The  above  derivation  is  dependent  only  on  the  assumption  of  quasi-equilibrium  for 
the  repelled  species  at  every  point  in  space. 

A.3.2  Current  Equation  for  Attracted  Species 

In  the  literature,  it  is  generally  assumed  that  the  plasma  species  that  is  attracted  to 
the  probe  experiences  a  "free-fall"  acceleration  toward  the  probe;  that  is,  the  kinetic 
energy  of  a  particle  at  a  position  where  the  electrostatic  potential  is  <Kx)  is  given  by 


m2,.  1  2  . ,  x 

j  v  (*)  =  2  mv.  +  q  <f>(x)  , 

where  v.  is  the  initial  velocity  of  the  particle  in  the  region  where  <Kx)  =  o  (Ref.  5). 


( A-10) 


For  an  initial  Maxwellian  distribution  of  particles  with  a  velocity  distribution  in  the  x 
direction  given  by 
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when  the  electrostatic  potential  is  zero,  the  average  velocity  of  all  the  particles  that  end 
up  moving  in  the  direction  of  the  acceleration,  when  the  local  <f>  is  non-zero,  is  given  by 
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The  first  term  in  Eq.  A-12  comes  from  the  particles  (half  the  total  initial  distribution) 
that  have  their  initial  vXq's  in  the  same  direction  as  their  final  vx's.  The  second  term 
comes  from  particles  whose  initial  vx0's  were  opposite  in  direction  to  the  final  vx's  but 
were  reversed  by  the  accelerating  fields. 

Similarly,  the  average  velocity  of  the  particles  whose  initial  vXq's  were  opposite  to 
the  accelerating  force  but  whose  x  velocities  were  never  reversed  in  direction  is  given  by 
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When  <{>  =  0, 
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(A-13) 


from  Eq.  A-4.  The  normalization  density  in  Eqs.  A-12  and  A-13  is  the  total  density  nQ. 
Thus,  the  flux  moving  in  either  the  positive  or  negative  direction  with  <j>  =  q  is  again 
(Hq/2)  vy^,  as  discussed  in  Section  A.3.1.  When  $  is  not  zero,  the  flux  in  the  direction  of 
the  accelerating  force  is  nQ  vx+ ,  and  in  the  direction  opposite  to  the  accelerating  force  is 

n0  vx_  • 

Equations  A-12  and  A-13  can  be  put  into  the  form 


where 


(B,v  ) 
0 


B  =  4|4'/kT  . 


( A-14) 


These  equations  can  be  numerically  integrated  very  simply  for  various  fixed  values  of  B. 
The  ratios  vx+  _/vj  are  plotted  vs.  B  in  Figure  25.  In  the  code,  these  functions  are 
approximated  numerically  to  give  the  flux  of  the  attracted  particles  out  of  each  mesh 
node  (with  potential  <)>)  in  the  direction  of  the  accelerating  force  (vx+)  and  opposite  the 
accelerating  force  (vx  ). 

It  is  interesting  that  vx^  can  be  approximated  fairly  well  by 
VT 

A  +  4  nq  4>/kT  . 

In  the  limit  of  large  4>,  this  velocity  approaches  /2q<J>/m;  that  is,  all  the  particles  are 
moving  with  the  velocity  corresponding  to  a  kinetic  energy  equal  to  the  potential  energy, 
Conversely,  at  large  values  of  <f>,  vx_  approaches  zero,  which  means  that  there  is  no 
flux  of  particles  opposing  the  accelerating  force. 
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Figure  25*  Average  particle  velocities  as  function  of  local  potential,  + 


AJJ  Plasma  Boundary  Conditions  at  Inner  Body  and  Outer  Conductor 

For  the  repelled  species,  the  flux  of  particles  into  the  inner  body  (that  is,  up  to  the 
body  surface,  where  they  are  assumed  to  be  captured  or  neutralized)  is  taken  to 
be  ng/2  vTxf  where  ng  is  the  total  density  of  particles  just  outside  the  body  (Ref.  5) 
and  Vy^  is  the  average  thermal  velocity  in  the  x  direction  of  half  the  particles  nB  (Eq. 
A-4).  This  boundary  condition  is  consistent  with  the  assumption  that  the  repelled  species 
is  everywhere  in  quasi-equilibrium;  that  is,  half  of  the  particles  are  moving  toward  the 
inner  body  and  half  are  moving  away,  even  adjacent  to  the  body.  In  reality,  right  at  the 
surface,  there  can  be  no  outward  current  due  to  shadowing  by  the  body  (in  the  absence  of 
secondary  emission,  which  has  been  ignored  thus  far  in  this  code). 
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However,  to  be  consistent  with  the  analytical  treatments  (Refs.  5,10)  and  the  quasi- 
equilibrium  assumption  for  the  repelled  particles,  shadowing  by  the  body  is  ignored. 

At  the  outer  boundary,  the  density  of  the  repelled  species  is  assumed  to  be  fixed  at 
the  bulk  value,  n^.  The  net  flux  of  these  particles  across  the  outer  boundary  is  determined 
by  Eq.  A-9,  where  n2  is  the  bulk  density  at  the  outer  boundary,  n^  is  the  density  at  the 
first  mesh  point  inside  the  boundary,  and  E  is  the  electric  field  between  these  two  points. 
The  magnitude  of  this  current  is  determined  self-consistently  in  the  code  to  satisfy  the 
time-dependent  rates  of  change  of  density  inside  the  volume  and  the  current  of  this 
species  into  the  center  body. 

For  the  attracted  species  the  flux  into  the  center  body  is  given  by  nRvx+/ 
where  ng  is  now  the  total  density  of  the  attracted  species  adjacent  to  the  center  body. 
Again,  everywhere  outside  the  center  body,  including  the  immediately  adjacent  mesh 
positions,  it  is  assumed  in  the  code  that  there  is  a  flux  of  particles  away  from  the  center 
body  given  by  nR  vx Very  close  to  the  center  body,  this  outward  flux  should  be  reduced 
by  shadowing.  However,  this  shadowing  effect  has  been  ignored  in  the  code.  Fortunately, 
for  large  attracting  potentials,  vx_  approaches  zero,  so  the  shadowing  effect  is  correctly 
accounted  for  at  large  attracting  potentials* 

At  the  outer  boundary,  the  density  of  the  attracted  species  is  also  fixed  at  the  value 
of  the  bulk  density,  nQ#  Thus,  since  the  potential  of  the  outer  boundary  is  taken  to  be 
zero,  there  is  an  inward  flux  of  attracted  particles  (nQ/2)vj^.  However,  there  is  also  an 
outward  flux  of  attracted  particles  from  the  first  mesh  point  inside  the  outer  boundary, 
given  by  n^  vx_ ,  where  n^  is  the  total  density  of  attracted  particles  at  the  first  zone 
inside  the  boundary.  Thus,  again,  the  net  flux  of  the  attracted  particles  across  the  outer 
boundary  is  determined  self-consistently  by  the  code  to  satisfy  the  time  rate  of  change  of 
the  particle  densities  in  the  volume  and  the  current  into  the  center  body. 

A3A  Simulation  of  Emission  Current 

At  present,  the  code  cannot  simulate  the  emission  of  a  discrete  beam  of  charge 
particles  with  finite  initial  energies.  The  simulation  of  emission  current  that  is  now  in  the 
code  is  equivalent  to  the  outward  emission  of  very  high-energy  particles,  approximately 
uniformly  from  the  surface  area  of  the  emitting  body.  The  high-energy  condition  means 
that  the  emitted  particles  are  slowed  down  negligibly  by  the  body  potential  and  their 
transit  time  from  the  emitting  body  to  the  outer  boundary  of  the  simulation  volume 
(essentially  a  few  Debye  lengths  away)  is  negligible  compared  to  the  response  times  of  the 


plasma*  This  condition  means  that  the  present  code  cannot  simulate  the  nonlinear  effects 
where  the  body  potential  would  approach  the  beam  energy  and  significantly  slow  down  the 
beam  particles* 

The  following  is  the  procedure  used  in  the  code.  At  the  start  of  a  plasma  problem, 
the  Poisson-solver  portion  of  the  code  is  run  first  for  a  unit  charge  differential  between 
the  emitting  body  and  the  outer  boundary.  The  cede  iterates  the  electrical  fields  until  V  x 
E  due  to  this  charge  differential  is  essentially  zero  throughout  the  volume,  and  then  stores 
these  normalized  fields  for  future  use*  In  the  plasma  problem,  the  amount  of  charge  AQ 
that  is  emitted  during  each  time  step  is  assumed  to  instantly  arrive  at  the  outer  boundary 
and  instantly  distribute  itself  over  the  outer  boundary,  consistent  with  a  quasi-static 
charge  differential.  The  electric  fields  at  every  point  inside  the  simulation  volume  are 
then  incremented  by  the  normalized  precharge  fields  calculated  at  the  beginning  of  the 
plasma  problem  times  the  charge  AQ  emitted  during  that  time  step.  The  total  resulting 
electric  fields  (the  previous  fields  plus  the  increments  due  to  AQ)  then  act  on  the  two 
plasma  species  during  that  time  step. 

A«4  CODE  STATUS  AND  PRELIMINARY  RESULTS 

As  mentioned  in  Section  A.1,  an  early  version  of  this  code,  which  treated  both  charge 
species  in  a  quasi-equilibrium  approximation,  was  used  to  try  to  reproduce  the  curve  in 
Reference  10  for  the  potential  of  a  non-emitting  probe  in  a  plasma  as  a  function  of  the 
ratio  of  the  electron  and  ion  temperatures.  In  Figure  26,  the  present  calculations  are 
compared  to  the  analytical  result  from  Reference  10.  Although  the  agreement  between 
the  two  results  is  outside  the  bounds  that  one  might  expect  just  due  to  numerical 
computational  techniques,  the  agreement  is  good  enough  to  indicate  considerable  promise 
for  the  code,  especially  since  the  two  models  do  not  use  precisely  the  same  physics  for  the 
attracted  species,  and  one  geometry  is  cylindrical  and  the  other  spherically  symmetric. 
Also,  it  should  be  pointed  out  that  at  least  part  of  the  discrepancy  at  the  larger  values  of 
T^Tj  results  from  these  calculations  being  terminated  before  the  slow-moving  ions  (small 
Tj)  had  reached  a  complete  steady-state  condition.  In  this  early  version  of  the  code, 
explicit  time  integrations  (using  the  velocities  at  the  beginning  of  the  time  steps)  were 
used  to  move  the  electrons  and  ions  during  the  time  steps.  When  the  electron  and  ion 
velocities  are  widely  different,  as  is  the  case  with  large  T^Tj,  computational  instabilities 
in  the  electron  densities  usually  result  if  the  time  step  is  made  larger  than  the  electron 
pt.iMitii  pevlfid  •«»  *pr»n«l  ii| i  Ihti  mh  iilrtlloiii  lot  I  lt«'  \li»wly  r ••%!>« unit ng  loin.  I  Hi  1 1  •  till  t«*% 


67 


such  as  these  can  be  overcome  with  implicit  integration  techniques  that  calculate  the 
change  in  carrier  densities  during  a  time  step,  using  some  average  velocity  during  that 
time  step.  This  type  of  integration  would  be  implemented  in  any  future  version  of  the 
code. 


Figure  26*  Comparison  of  present  code  with  analytical  predictions  from 
Reference  10  (no  secondary  emission  from  probe).  Te  and  Tj 
are  temperatures  of  electrons  and  ions,  respectively. 


The  present  status  of  the  code  is  that  the  free-fall  approximation  for  the  attracted 
species  has  been  programmed  in  a  particular  manner  and  a  crude  form  of  implicit 
integration  has  been  implemented.  However,  there  are  some  computational  and  physics 
difficulties  with  this  formulation  that  require  some  additional  work.  It  is  anticipated  that 
a  more  stable  implicit  integration  procedure  and  the  required  new  physics  formulation 
could  be  accomplished  with  a  modest  effort. 
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A*5  FUTURE  CODE  CAP ABL ITIES 


There  are  some  additional  problems  for  which  a  code  such  as  the  present  one  would  be 
useful,  but  they  would  require  some  additional  modifications  to  the  present  coding.  Two 
of  the  more  important  applications  are  avalanche  ionization  of  the  neutral  background  gas 
by  energetic  electrons,  either  those  emitted  from  the  probe  or  those  returning,  and  the 
effects  of  photo- ionization  and  the  resulting  photon-driven  current  when  a  photon  beam 
interacts  with  the  background  neutral  gas. 

In  both  these  problems,  the  question  is  how  to  combine  the  newly  created  electrons 
and  ions,  which  would  have  some  velocity  distributions,  with  the  already  existing  electron 
and  ion  densities,  which  would  normally  have  different  distributions. 

A  simple,  but  crude,  approach  would  be  to  assume  that  the  newly  created  particles 
have  the  same  velocity  distributions  as  the  previously  existing  particles  at  every  location 
in  space.  This  approach  might  not  be  too  erroneous  for  avalanche  ionization,  where  the 
dominant  velocities  for  the  new  and  old  electrons  would  at  least  be  pointing  in  the  same 
directions.  However,  it  would  be  very  poor  for  trying  to  simulate  a  photon-driven  current 
at  some  angle  to  the  dominant  velocities. 

A  somewhat  more  refined  approach  would  be  to  combine  the  momenta  of  the  new  and 
old  particles  vectorial ly  at  each  mesh  location  and  use  the  combined  momentum 
(magnitude  and  direction)  to  specify  the  velocity  of  the  shifted  Maxwellian.  Obviously, 
this  approach  is  not  vigorously  correct  either.  However,  Crevier  (Ref.  22)  has  shown  that 
a  similar  procedure,  when  used  with  SGEMP  electrons,  gave  results  that  agreed  very  well 
with  equivalent  particle-pushing  techniques.  Apparently  the  details  of  the  distributions 
are  not  too  critical,  at  least  for  gross  effects,  as  long  as  the  momentum  and  energy  of  the 
plasma  clouds  are  conserved.  However,  adoption  of  this  approach  would  require  a  mod¬ 
ification  of  the  free-fall  formalism,  since  the  free-fall  approach  assumes  that  the  net 
velocity  of  the  plasma  cloud  is  a  function  only  of  the  local  potential.  A  form  of  the  free- 
fall  coding  which  has  been  developed  recently  would  be  compatible  with  this  approach  for 
including  electron  and  photon  ionization.  In  addition,  it  appears  that  the  new  coding 
method  wilt  be  computationally  more  stable.  Again,  these  modifications  to  the  code 
should  not  require  an  extensive  effort. 
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